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Three  Dimensional  Transient  Analysis  of  Microstrip  Circuits  in  Multilayered  Anisotropic  Media 


Under  the  sponsorship  of  the  ONR  contract  N00014-90-J-1002  we  have  published  9 
refereed  journal  and  conference  papers. 

The  coupled-wave  theory  is  generalized  to  analyze  the  diffraction  of  waves  by  chiral 
gratings  for  arbitrary  angles  of  incidence  and  polairizations.  Numerical  results  for  the 
Stokes  parameters  of  diffracted  Floquet  modes  versus  the  thickness  of  chiral  gratings  with 
various  chiralities  are  calculated.  Both  horizontal  and  vertical  incidences  are  considered  for 
illustration.  The  diffracted  waves  from  chiral  gratings  eure  in  general  elliptically  polarized; 
and  in  some  p<irticular  instances,  it  is  possible  for  chiral  gratings  to  convert  a  linearly 
polarized  incident  field  into  two  nearly  circularly  polarized  Floquet  modes  propagating  in 
different  directions. 

A  general  spectral  domain  formulation  to  the  problem  of  radiation  of  arbitruy 
distribution  of  sources  embedded  in  a  horizontally  stratified  arbitrary  magnetized  linear 
plasma  is  developed.  The  fields  are  obtained  in  terms  of  electric  and  magnetic  type  dyadic 
Green’s  functions.  The  formulation  is  considerably  simplified  by  using  the  kDB  system 
of  coordinates  in  conjunction  with  the  Fourier  transform.  The  distributional  singul£ir 
behavior  of  the  various  dyadic  Green’s  functions  in  the  source  region  is  investigated  ^Lnd 
taken  into  account  by  extracting  the  delta  function  singiilarities.  Finally,  the  fields  in  any 
arbitrary  layer  are  obtained  in  terms  of  appropriately  defined  global  upward  and  downward 
reflection  and  treinsmission  matrices.  ^ 

We  investigated  a  method  for  the  calculation  of  the  current  distribution,  resistance, 
and  inducteince  matrices  for  a  system  of  coupled  superconducting  transmission  lines  having 
finite  rectangular  cross  section.  These  calculations  allow  accurate  characterization  of  both 
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high-Tt  and  low-Te  superconducting  strip  transmission  lines.  For  a  single  stripline  geometry 
with  finite  ground  planes,  the  current  distribution,  resistance,  inductance,  and  kinetic 
inductance  are  calculated  as  a  function  of  the  penetration  depth  for  various  film  thickness. 
These  calculations  ue  then  used  to  determine  the  penetration  depth  for  Nb,NbN,  and 
YBa2CuaOj^a  superconducting  thin  films  from  the  measured  temperature  dependence  of 
the  resonant  frequency  of  a  stripline  resonator.  The  calculations  are  also  used  to  convert 
measured  temperature  dependence  of  the  quality  factor  to  the  intrinsic  surface  resistance 
as  a  function  of  temperature  for  a  Nb  stripline  resonator. 

The  electromagnetic  radiation  from  a  VLSI  chip  package  and  heatsink  structure 
is  analysed  by  means  of  the  finite-difference  time-domain  (FD-TD)  method.  The  FD- 
TD  algorithm  implemented  incorporates  a  multi-zone  gridding  scheme  to  accommodate 
fine  grid  cells  in  the  vicinity  of  the  heatsink  and  package  cavity  and  sparse  gridding  in 
the  remainder  of  the  computational  domain.  The  issues  pertaining  to  the  effects  of  the 
heatsink  in  influencing  the  overall  radiating  capacity  of  the  configuration  are  addressed. 
Analyses  are  facilitated  by  using  simplified  heatsink  models  and  by  using  dipole  elements  as 
sources  of  electromagnetic  energy  to  model  the  VLSI  chip.  The  potential  for  enhancement 
of  spurious  emissions  by  the  heatsink  structure  is  examined.  For  heatsinks  of  typical 
dimensions,  resonance  is  possible  within  the  low  gigahertz  frequency  range. 

Because  the  effects  of  diffraction  during  proximity-print  x-ray  lithography  are  of 
critical  importance,  a  number  of  previous  researchers  have  attempted  to  calculate  the 
diffraction  patterns  and  minimum  achievable  feature  sizes  as  a  function  of  wavelength  and 
gap.  Work  to  date  has  assumed  that  scalar  diffraction  theory  is  applicable-as  calculated,  for 
example,  by  the  Rayleigh-Sommerfeldformulation-andthat  Kirchhoff  boundary  conditions 
can  be  applied.  Kirchhoff  boundary  conditions  assume  that  the  fields  (amplitude  and 
phase)  are  constant  In  the  open  regions  between  absorbers,  and  a  different  constant  in 
regions  just  under  the  absorbers  (i.e.,  that  there  are  no  fringing  fields).  An  x-ray  absorber 
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is,  however,  best  described  as  a  lossy  dielectric  that  is  tens  or  hundreds  of  wavelengths 
tall,  and  hence  Kirchhoff  boundary  conditions  are  unsuitable.  We  have  instead  used  two 
numerical  techniques  to  calculate  accurate  diffracted  fields  from  gold  absorbers  for  two 
ceises:  a  30  nm-wide  line  at  A  =  4.5  nm,  and  a  100  nm-wide  line  at  A  =  1.3  nm.  We 
show  that  the  use  of  Kirchhoff  bound2Lry  conditions  introduces  unphysically  high  spatial 
frequencies  into  the  diffracted  fields.  The  suppression  of  these  frequencies-which  occurs 
naturally  without  the  need  to  introduce  an  extended  source  or  broad  spectrum-improves 
exposure  latitude  for  mask  features  near  0.1  ftm  and  below. 

In  order  to  understand  the  physical  meaning  of  rational  reflection  coefficients  in 
one-dimensional  inverse  scattering  theory  for  optical  waveguide  design,  we  have  studied 
the  relation  between  the  poles  of  the  transverse  reflection  coefficient  and  the  modes  in 
inhomogeneous  dielectrics.  By  using  a  stratified  medium  model  it  is  shown  that  these 
poles  of  the  reflection  coefficient  have  a  one-to-one  correspondence  to  the  discrete  modes, 
which  are  the  guided  and  leaky  modes.  The  radiation  modes  have  continuous  real  values  of 
transverse  wave  numbers  and  are  not  represented  by  the  poles  of  the  reflection  coefficient. 
Based  on  these  results,  applications  of  the  Gel’fand-Levitan-Marchenko  theory  to  optical 
waveguide  synthesis  with  the  rational  function  representation  of  the  transverse  reflection 
coefficient  are  investigated. 

In  compact  modules  of  high  performance  computers,  signal  transmission  lines  be¬ 
tween  integrated  circuit  chips  are  embedded  in  multilayered  dielectric  medium.  These 
signal  lines  are  usually  placed  in  different  layers  and  run  perpendicular  to  each  other.  The 
interaction  between  the  orthogonal  crossing  lines  and  the  signal  line  affects  its  propagation 
characteristics  and  may  cause  considerable  signal  distortion. 
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The  interaction  of  a  pair  of  crossing  lines  in  isotropic  medium  has  been  studied 
using  a  time-domjun  approach,  where  coupling  is  described  qualitatively.  This  method 
becomes  computationally  expensive  when  the  number  of  crossing  lines  increases.  With 
many  identical  crossing  strips  uniformly  distributed  above  the  signal  line,  the  transmis¬ 
sion  properties  are  chzu'acterized  by  stopbands  due  to  the  periodicity  of  the  structure. 
Periodic  structure  have  been  investigated  using  frequency-domain  methods.  Periodically 
nonuniform  microstrip  lines  in  an  enclosure  are  analyzed  on  the  b^lsis  of  a  numerical  field 
calculation.  A  technique  based  on  the  network-analytical  formulism  of  electromagnetic 
fields  has  been  used  to  analyze  striplines  and  finlines  with  periodic  stubs.  The  propagation 
cheiracteristics  of  waves  along  a  periodic  array  of  parallel  signal  lines  in  a  multilayered 
isotropic  structure  in  the  presence  of  a  periodically  perforated  ground  plame  and  that  in 
a  mesh-plane  environment  have  been  studied.  More  recently,  the  effect  of  the  geometrical 
properties  on  the  propagation  characteristics  of  strip  lines  with  periodic  crossing  strips 
embedded  in  a  shielded  one-layer  isotropic  medium  have  been  investigated.  Both  open 
and  closed  multilayered  uniaxially  anisotropic  structures  are  considered.  A  full-wave  anal¬ 
ysis  is  used  to  study  the  propagation  characteristics  of  a  microstrip  line  in  the  presence 
of  crossing  strips.  The  signal  line  and  the  crossing  strips  are  assumed  to  be  located  in 
two  arbitrary  layers  of  a  stratified  uniaxially  anisotropic  medium.  An  integral  equation 
formidation  using  dyadic  Green’s  functions  in  the  periodicsdly  loaded  structure  is  derived. 
Galerkin’s  method  is  then  used  to  obtain  the  eigenvalue  equation  for  the  propagation  con¬ 
stant.  The  effects  of  anisotropy  on  the  stopband  properties  are  investigated.  Numerical 
results  for  open  and  shielded  three-layer  uniaxially  anisotropic  media  are  presented. 

For  microwave  integrated  circuit  applications,  the  characteristics  of  interconnects 
have  been  investigated  for  the  propagation  modes,  time  response,  crosstalk,  coupling, 
delay,  etc.  In  these  analyses,  it  is  assumed  that  quasi-TEM  modes  are  guided  along  the 
multiconductor  transmission  lines.  The  analysis  were  performed  for  arbitrary  number  of 
transmission  lines  where  the  load  and  the  source  conditions  were  presented  in  terms  of  the 
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modal  reflection  and  transmission  coeflicient  matrices. 


To  perform  the  quasi-TEM  analysis,  the  capacitance  matrix  for  the  miilticonductor 
transmission  line  has  to  be  obtained  first.  Both  the  spectral  and  the  spatial  domain 
methods  have  been  proposed  to  calcvilate  the  capacitance  matrix.  In  the  spectral  domain 
methods,  two  side  walls  are  used  to  enclose  the  whole  transmission  line  structure,  and  the 
thickness  of  the  strip  lines  has  not  been  considered.  In  using  the  spatial  domain  method, 
the  structure  has  to  be  truncated  to  a  finite  extent  to  medce  the  numerical  implementation 
fezisible.  The  infinite  extent  of  the  structure  was  also  incorporated,  but  only  a  two-layer 
medium  was  considered. 

In  practiced  microwave  integrated  circuits,  the  dielectric  loss  due  to  the  substrate 
and  the  conductor  loss  due  to  the  metallic  strips  are  also  studied  in  the  einalysis  of  circuit 
performances. 

Based  on  the  scalar  Green’s  ftmction,  a  set  of  coupled  integral  equations  is  obtained 
for  the  charge  distribution  on  the  strip  surfaces.  Pulse  basis  functions  and  a  point-matching 
scheme  is  used  to  solve  numerically  the  set  of  integral  equations  for  the  charge  distribution, 
and  hence  the  capacitance  matrix.  The  duality  between  the  electrostatic  formulation  amd 
the  magnetostatic  one  is  applied  to  calculate  the  inductance  matrix.  The  conductance 
matrix  is  obtained  by  using  the  duality  between  the  electrostatic  problem  and  the  current 
field  problem.  A  perturbation  method  is  used  to  calculate  the  resistance  matrix. 

Finally,  a  transmission  line  analysis  is  derived  to  obtain  the  transfer  matrix  for  multi- 
conductor  uniform  lines,  which  significantly  reduces  the  effort  in  treating  the  load  and  the 
source  conditions.  Transient  responses  are  obtained  by  using  the  Fourier  transform.  The 
results  for  two  coupled  lines  are  obtained. 

With  the  ever  increasing  speed  and  density  of  modern  integrated  circuits,  the  need 
for  electromagnetic  wave  analysis  of  phenomena  such  as  the  propagation  of  transient  sig- 
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nals,  especially  the  distortion  of  signal  pulses,  becomes  crucial.  One  of  the  most  important 
causes  of  pulse  distortion  is  the  frequency  dependence  of  conductor  loss,  which  is  caused 
by  the  “skin  effect”,  and  which  can  be  incorporated  into  the  circuit  models  for  transmis¬ 
sion  lines  as  frequency-dependent  resistance  and  inductance  per  unit  length.  Efficient  and 
accurate  algorithms  for  calculating  these  parameters  are  increasingly  important. 

We  have  developed  a  hybrid  cross-section  finite  element /coupled  integral  equation 
method.  The  technique  is  a  combination  of  a  cross-section  finite  element  method,  which  is 
best  for  high  frequencies.  An  interpolation  between  the  results  of  these  two  methods  gives 
very  good  results  over  the  entire  frequency  range,  even  when  few  basis  functions  are  used. 

In  the  cross-section  method,  we  divide  each  conductor  into  triangular  patches  and 
choose  one  of  the  patches  from  the  return  conductor  to  be  our  reference.  We  then  calculate 
the  resistance  and  inductance  matrices  for  the  patches.  Using  two  conditions  on  the  system, 
that  the  total  current  in  each  wire  is  the  sum  of  the  currents  in  the  patches,  and  that  the 
voltage  on  each  patch  in  a  wire  must  be  the  same  (no  transverse  currents),  we  can  reduce 
the  matrices  for  the  patches  to  the  matrices  for  the  wires.  In  the  Weeks  method,  the 
patches  are  rectangles,  and  the  quadruple  integral  is  done  quite  easily  in  closed  form. 
However,  it  is  also  possible  to  evaluate  the  quadruple  integral  in  closed  form  for  triangular 
patches,  although  the  mathematics  leading  to  this  result  is  quite  involved,  and  the  final 
form  of  the  answer  is  complicated.  We  therefore  use  triangular  patches  as  the  most  flexible 
means  of  modelling  conductors  with  arbitrary  cross-sections;  polygons  are  covered  exactly, 
and  we  are  able  to  model  quite  closely  other  shapes,  such  as  circles. 

As  frequency  increases,  the  need  to  keep  the  uniform  current  approximation  valid  in 
the  patches  requires  either  the  addition  of  many  more  patches  as  the  skin  depth  decreases, 
or  a  redistribution  of  the  existing  patches  to  the  surface,  where  the  current  is.  However, 
changing  the  distribution  of  patches  makes  it  necessary  to  recalculate  the  resistance  and 
inductance  matrices  of  the  patches,  thus  increasing  the  computation  time.  Since  we  use  a 
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surface  integral  equation  method  for  high  frequencies,  we  do  not  change  the  distribution 
of  the  triangular  patches  for  the  cross-section  method  as  we  increase  the  frequency. 

For  high  frequencies,  we  use  a  coupled  surface  integral  equation  technique.  Under 
the  queisi-TEM  assumption,  the  frequency-dependent  resistance  and  inductance  result  from 
the  power  dissipation  and  magnetic  stored  energy,  which  can  be  calculated  by  solving  a 
magnetoquasistatic  problem,  with  the  vector  potential  satisfying  Laplace’s  equation  in  the 
region  outside  all  the  conductors.  The  resistance  and  inductance  are  usually  given  by 
integrals  of  these  field  quemtities  over  the  cross-sections  of  the  wires,  but  by  using  some 
vector  identities  it  is  possible  to  convert  these  expressions  to  integrals  only  over  the  surfaces 
of  the  wires.  These  expressions  contain  only  the  current  at  the  surface  of  each  conductor, 
the  derivative  of  that  current  normal  to  the  surface,  and  constants  of  the  vector  potential. 
A  coupled  integral  equation  is  then  derived  to  relate  these  quantities  through  Laplace’s 
equation  and  its  Green’s  function  outside  the  conductors  and  the  diffusion  equation  and  its 
Green’s  function  inside  the  conductors.  The  method  of  moments  with  pulse  basis  functions 
is  used  to  solve  the  integral  equations.  This  method  differs  from  previous  work  in  that  the 
calculation  of  resistance  and  inductance  is  based  on  power  dissipation  and  stored  magnetic 
energy,  rather  than  on  impedance  ratios.  It  will  therefore  be  more  easily  extended  to 
structures  where  non-TEM  propagation  can  occur. 

For  the  intermediate  frequency  range,  where  the  conductors  are  on  the  order  of  the 
skin  depth,  were  foimd  it  very  efficient  to  interpolate  between  the  results  of  the  cross- 
section  and  surface  methods.  The  interpolation  function  was  based  on  the  average  size  of 
the  conductors,  measured  in  skin  depths,  and  was  of  the  form  1/(1  -|-  0.16o’/5*),  where  it  a 
is  the  average  cross-section  of  the  conductors,  and  6  is  the  skin  depth. 
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The  Propagation  Characteristics  of  Signal  Lines  with  Crossing 
Strips  in  Multilayered  Anisotropic  Media 

C.  W.  Lam,  S.  M.  Ali,  and  J.  A.  Kong 

Department  of  Electrical  Engineering  and  Computet  Science 
and  Research  Laboratory  of  Electronics 
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Abstract —  In  this  paper,  full  modal  analysis  is  used  to  study  the  dispersion  character¬ 
istics  of  microstrip  lines  periodically  loaded  with  crossing  strips  in  a  stratified  uniaxiaUy 
anisotropic  medium.  Dyadic  Green's  functions  in  the  spectral  domain  for  the  multilayered 
medium  in  conjunction  with  the  vector  Fourier  transform  (VFT)  are  used  to  formulate  a 
coupled  set  of  vector  integrsJ  equations  for  the  current  distribution  on  the  signal  line  and 
the  crossing  strips.  Galerkin’s  procedure  is  applied  to  derive  the  eigenvalue  equation  for 
the  propagation  constant.  The  effect  of  anisotropy  for  both  open  and  shielded  structures 
on  the  stopband  properties  is  investigated. 

I.  INTRODUCTION 

In  compact  modules  of  high  performance  computers,  signal  transmission  lines 
between  integrated  circuit  chips  are  embedded  in  multilayered  dielectric  media. 
These  signal  lines  are  usually  placed  in  different  layers  and  run  perpendicular 
to  each  other.  The  interaction  between  the  orthogonal  crossing  lines  and  the 
signal  line  affects  its  propagation  characteristics  and  may  cause  considerable  signal 
distortion. 

The  interaction  of  a  pair  of  crossing  lines  in  an  isotropic  medium  hr.s  been  stud¬ 
ied  using  a  time-domain  approach  [1],  where  coupling  is  described  qualitatively. 
This  method  becomes  computationally  expensive  when  the  number  of  crossing 
lines  increases.  With  many  identical  crossing  strips  uniformly  distributed  above 
the  signal  line,  the  transmission  properties  are  characterized  by  stopbands  due  to 
the  periodicity  of  the  structure.  Periodic  structures  have  been  investigated  using 
frequency-domain  methods.  In  [2],  periodically  nonuniform  microstrip  lines  in  an 
enclosure  are  analyzed  on  the  basis  of  a  numerical  Held  calculation.  A  technique 
based  on  the  network-analytical  formulism  of  electromagnetic  fields  has  been  used 
to  analyze  stripb'nes  and  finlines  with  periodic  stubs  [3].  The  propagation  charac¬ 
teristics  of  waves  along  a  periodic  array  of  parallel  signal  lines  in  a  multilayered 
isotropic  structure  in  the  presence  of  a  periodically  perforated  ground  plane  is 
studied  in  [4]  and  that  in  a  mesh-plane  environment  is  studied  in  [5].  More  re¬ 
cently,  the  effect  of  the  geometrical  properties  on  the  propagation  characteristics  of 
strip  lines  with  periodic  crossing  strips  embedded  in  a  shielded  one-layer  isotropic 
medium  have  been  investigated  (6). 
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In  this  paper,  both  open  and  closed  multilayered  uniaxially  anisotropic  struc¬ 
tures  are  considered.  A  full-wave  analysis  is  used  to  study  the  propagation  char¬ 
acteristics  of  a  microstrip  line  in  the  presence  of  crossing  strips.  The  signal  line 
and  the  crossing  strips  are  assumed  to  be  located  in  two  arbitrary  layers  of  a 
stratified  uniaxially  anisotropic  medium.  An  integral  equation  formulation  using 
dyadic  Green’s  functions  in  the  periodically  loaded  structure  is  derived.  Galerkin’s 
method  is  then  used  to  obtain  the  eigenvalue  equation  for  the  propagation  con¬ 
stant.  The  effects  of  anisotropy  on  the  stopband  properties  are  investigated. 
Numerical  results  for  open  and  shielded  three-layer  uniaxially  anisotropic  media 
asc  presented. 

n.  FORMULATION  OF  THE  PROBLEM 

In  this  section,  we  present  a  dyadic  Green’s  function  formulation  of  the  problem 
shown  in  Fig.  1(a)  where  the  microstrip  line  and  the  crossing  strips  are  placed 
at  two  different  interfaces  of  a  uniaxially  anisotropic  multilayered  medium.  The 
crossing  strips  are  assumed  to  be  placed  in  a  layer  (i)  and  the  signal  line  to  be 
in  a  layer  (j).  The  crossing  strips  are  considered  to  be  periodic  with  period  p  as 
shown  in  Fig.  1(b).  In  general,  the  permittivity  and  permeability  tensors  of  an 
arbitrary  layer  (/)  are  assumed  to  be  given  by 


( 

0 

^  \ 

h  = 

(1) 

Vo 

0 

«lz/ 

(  H 

0 

(2) 

%  =  1 

^  0 

0 

where  I  =  0, 1,2, ...,n, ...,t. 

For  the  stratified  medium,  the  electric  fields  in  layers  (»)  and  (j)  due  to  current 
distributions  Ji(r*)  and  Jj{^)  may  be  expressed  as 

£i(r)  =  xu;  dV  Wii{f,T’) .  7i(r')  iu;  J j dV  (3a) 

£y(f )  =  dV  Wjiif,  f')  .  7i(r')  +  ‘<^  ///^  dV  f’)  .  (3a) 

where  Gi^{r,f*)  is  the  dyadic  Green’s  function  in  layer  (f)  due  to  current  sources 
in  layer  (m). 

For  the  multilayered  structure  shown  in  Fig.  1(a),  the  current  distributions  on 
the  conducting  strips  are  assumed  to  be  surface  currents  lying  on  planes  transverse 
to  z.  Thus,  if  we  are  interested  in  the  transverse  electric  fields  and  in 
layers  (»)  and  (j),  respectively,  we  can  write 


(4a) 

(45) 
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Figure  1.  Geometrical  configuration  of  a  signal  line  in  layer  j  loaded  with 
crossing  strips  in  layer  t. 

where  is  the  transverse  electric  field  in  layer  (f)  due  to  current  sources  in 
layer  (m),  and  is  given  by 

/ fj^  .  KmiT')  (5) 

_  =(,)  _ 

Km{f)  18  the  surface  current  distribution  in  layer  (m),  and  Gi^{r,T')  is  the  (2x2) 

transverse  part  of  the  dyadic  Green’s  function 

Since  the  structure  is  assumed  to  be  periodic  in  the  y-direction,  the  electric  field 
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'E\*^  and  the  surface  currents  Km  can  be  expressed  using  the  Floquet  harmonic 
representation  in  the  y-direction.  In  this  case  we  can  get 

^ImC)  = 

where  I  =  t,j,  and 


Fik,n,rs)  =  ^ 

*jn 


kx  0n 
0n  ~kx 


K  F{k,a,  -K)  •  ^jiK) 

The  matrix  iij{k$n)  is  given  by 


yTM 


xTE 

'»y  j 


r  dkx  F(fc.„,f.)  .  .  Ki{k,n) 

B=-00’'-«> 


B=— OO 


oo  _ 


r^'i 

'p/2  ^g'\ 

dy'^rl{f,f>).Kmif') 

(6a)  - 

J— oo  •/- 

-pl2 

«>  fOO 

t 

(66) 

r  -a  a  o 

r,  =  ix  +  yy,  *,n  =  *fcx  +  yPn,  Pn  =  Po  + - 

P 

Here  p  is  the  period,  0q  is  the  propagation  constant  of  the  donoinant  harmonic 
in  the  Floquet  representation,  and  ^l^^^***’*’**)  spectral  dyadic  Green’s 

function.  _ ^  ^  ^  ^  ^  ^ 

Using  the  explicit  expressions  for  the  dyadic  Green’s  functions  G,-,-  ,  Gjy  ,  Gjj 

and  [7],  the  transverse  electric  fields  on  the  surface  of  the  conducting  strips 
in  layer  (i)  due  to  the  currents  in  layer  {j)  can  be  expressed  in  the  following  form 

AOO 

=  E  /  df‘^^(^»n.u)-Jij(hn)-Kj(k,n)  (7) 

P  n=-oo''-‘^ 

where  F{ktn,^$)  is  the  kernel  of  the  vector  Fourier  transform  (VFT)  given  by  [8] 


(8) 


and  Kj{kgn)  is  the  vector  Fourier  transform  of  the  surface  current  Kj(f'g).  It  is 
given  by 

-L  1  r  roo  =  _ 

(9) 


(10) 


whose  elements  for  different  t  and  j  are  given  in  Appendix  A. 

In  the  above,  the  transverse  electric  field  expressions  (I,m  =  i,j)  satisfy 
the  boundary  conditions  at  the  dielectric  interfaces  of  the  layered  medium.  Apply¬ 
ing  the  final  boundary  condition  that  the  tangential  electric  field  vanishes  on  the 
conducting  strips,  we  can  get  the  following  set  of  dual  vector  integral  equations 
for  the  currents  on  the  metallic  strips 


dkx  F{ktn,r$)  '(ij{ktn)  •Kj{kgn)  —  0,  Vj  G  Si  (H) 
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on  the  crossing  strips,  and 

oo  /OO 

^  ]  I  F{ktfi,rf)  •  ^j^[kgn)  ’  Ki{kg'n) 

n=-oo''-«» 

oo  OQ 

+  ^  ]  I  <ikx  F{kgn,rg)  •  ^jj{^kgn)  •  Kj{kgn)  =  0,  TgESj  (12) 
n=-oo''-'» 
on  the  signal  line. 

The  next  step  is  to  solve  this  coupled  set  of  vector  integral  equations  to  find 
the  dispersion  relation  for  the  signal  line  in  the  presence  of  the  crossing  strips. 

m.  GALERKIN’S  METHOD  AND  THE  EIGENVALUE  EQUATION 

The  formulation  up  to  this  stage  is  exact.  We  now  solve  the  set  of  vector  integral 
equations  (11)  and  (12)  by  using  Galerkin’s  method.  The  unknown  current  distri¬ 
butions  on  the  crossing  strips  ^^^(r,)  and  on  the  signal  line  Kj(rg)  are  expanded 
in  terms  of  the  appropriate  vector  basis  functions  as  follows: 

_  M  R  _  _ 

^mT(®,y)  •  j4-mT  (13) 

»n=l  r=l 

_  K  Qi  ^ 

E  (14) 

fc=l9=-Ql 

where  Ki{x,y)  and  Kj(x,y)  are  the  surface  currents  on  the  crossing  strips  and 

^e  signal  line,  respectively,  $mr  and  are  the  basis  functions,  Amr  and 

B/ig  are  the  expansion  coefficients. 

U^g  (13)  and  (14),  the  vector  Fourier  transform  (VFT)  of  the  currents  Ki{rg) 
and  Kj{fg)  are  obtained  as 

M  R  _ 

^iikgn)~  E  E  ^l+)mT.n{kxt0)  '  Amr  (15) 

m=l r=l 

^  _  K  Qi  ^  _ 

A^,(fc*n)  =  E  E  ^[+)hq,n{kx,0)'Bt^  (16) 

k=lq=-Qy 

where 

~  (2n)^  JJ  •  ^inr(*,y)  (17) 

and 

=  (2^  ^(**’*’  (18) 
Substituting  (15)  and  (16)  into  (11)  and  (12),  we  obtain 

y  '  y  1  /  ^kx  F{kgn,rg)  •  ^^^[kgn)  •  U (kx,0)  •  Amr 

m=lr=ln=-oo''~®° 
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K  Q2  oo 


fc=l  5=:-<Jl  n=-oo  •' 


=  0 


(19) 


for  r,  e  5i,  and 
M  R  00 


EE  E  /  dfc*  ■F(kttn’‘t)  ‘  (ji(^sn)  • 

m=lr=ln=— 00*'“®°  ’ 

K  Ql  00  -00  _ _  =  _  =  _ 

fc=l  9=-Qi  n=-oo  •'"** 


=  0 


(20) 


for  r,  6  5j. 

— ^  - 

Multiplying  (19)  by  ^u»(*,  y)  and  integrating  over  the  support  of  A’,(f,)  for 

u  =  1,2,. ..,M  and  v  =  we  obtain 

m=lr=ln=-oo*'-«> 
t=l5=-Ql  ns-oo*'  ®° 

(21) 

Similarly,  multiplying  (20)  by  ’J',(z)e~*^y  and  integrating  over  the  support  of 
Kj(ft)  for  s  =  1,2,, ..,Jf  and  i  =  —<?!).•. ,0,-*.,Q2>  obtain 

EE  E 

m=lr=ln=-oo''-®® 

*=l,=-<3in=-oo-'-«> 

(22) 

Equations  (21)  and  (22)  constitute  a  system  of  (5  +  T)  linear  algebraic  equations 
with  5  =  MR  and  T  =  K(Qi  +  Q2  +  f ))  and  may  be  written  in  matrix  form  as 

^•c  =  0 

where 

l^lllsxS  [•^12]sxr 
.(■^2l]rxS  [•^22jrxTj 


N  = 


(23o) 

(236) 


and 


MmrJsxl 

'  ~  [-PfcalTxl  J 


(23c) 


The  Propagation  Charaeterietiee  of  Signal  Linet  with  Croiting  Strips 


1011 


I 

k  oc  C< 
f  cmcc 
K  oc  Cl 

c  •  c:  c 

L  OC  Cl 
i  C  •  C  C 

ic  OC  C  I 


c  o  c  c 

a,  < 

L  Oi  C 

Cl.  i 

(  Oi 
Cl  i 

‘  O-  ' 

,  C  4.  w  < 


Each  element  of  the  submatrices  of  N  is  given  by 

po  #oo  ^ 

n=-oo''-®° 

(24a) 

~  X]  /  ^(-)uv  ■  ^(+)jbo 

n=-oo*'-<»  ’  *’ 

(24b) 

PQ  OP 

n=-oo’'-‘» 

(24c) 

po  OP 

l^22W,=  E  / 

n=-oo*'-°°  ’ 

(24(f) 

For  nontrivial  solution  to  exist,  the  determinant  of  the  coefficient  matrix  of 
(21)  and  (22)  must  be  zero, 

det[^(a.,^)]  =0  (25) 

This  is  the  eigenvalue  equation  for  the  propagation  constant  P  which  describes 
the  dispersion  relation  of  the  loaded  microstrip  line  in  the  multilayered  anisotropic 
medium. 

_ The  next  step  is  to  choose  appropriate  basis  functions  for  the  surface  currents 

Kj{fi)  and  Ki{rt)  on  the  signal  line  and  the  crossing  strips,  respectively.  The 
expansion  functions  we  use  are 

«mr(*,Vj-2^[  0  (?m(x,I<)J  [  0  i?r(y,u;0j 

and 

^*^*^"27r[  0  gfc(x,ta,)J  (27) 


where 


1  .  2n7rQ 

•fn(Q:,7)  =  -  sin - 

7  7 


(J.(<.,7)  =  -  (29) 

(30) 

7  7 

Wi  and  wj  are  the  widths  of  the  crossing  strips  and  the  signal  line,  respectively, 
and  Li  is  the  length  of  the  crossing  strips.  When  choosing  the  basis  functions 
for  the  surface  currents,  it  should  be  borne  in  mind  that  the  current  cannot  have 
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a  normal  component  to  the  strip  edges.  Futhermore,  the  edge  condition  for  the 
parallel  component  should  be  considered.  By  substituting  (26)  and  (27)  into  (17) 
and  (18),  respectively,  we  can  get 


and 


y^lmr.n 


M) 


V 

(2ff)3fc,„ 


±Jfm(fc*,I<)  0  1 

0  rm(fc*,ii)J 

lr(^n,tft)  0 

®  ^(±),(^»»,W«) 


P^tpi  ±Xit(kx,Wj)  0 

(2ir)3jb,„  0  Yk{kx,wj) 


(31) 


(32) 


where 


Yn{an)  =  |{  JoI(n  -  l)7r  +  ^]  +  Jo[(n  -  1)^  -  ^]} 

_  ,  .  — t  r  ,  niT  cty . 

-  0272  “7)«n(—  T  — ) 


—  c 


■(nirTa7)Mn(Y=^^)] 


(33) 

(34) 

(35) 


Equations  (31)  and  (32)  are  then  substituted  into  the  determinantal  equation  (25) 
for  the  calculation  of  the  dispersion  characteristics. 


IV.  NUMERICAL  TREATMENT  AND  RESULTS 

In  this  section,  we  present  numerical  results  for  open  and  closed  three-layer  struc¬ 
tures  with  the  crossing  strips  and  the  signal  line  embedded  in  two  different  layers 
as  shown  in  Figs.  2(a)  and  (b).  In  numerical  calculation,  the  infinite  series  of 
Floquet  modes  and  the  basis  functions  are  truncated.  The  ranges  of  indices  in 
(24)  are  chosen  as;  —10  <n<9,  fc  =  l,  — 1<9<0,  l<m<3,  r  =  l,s  = 
1,  -1  <  t  <  0,  1  <  u  <  3,  and  r  =  1.  It  can  be  seen  that  each  element  in  the  co¬ 
efficient  matrix  can  be  reduced  to  a  sum  of  TE  and  TM  terms,  a  summation  over 
n  Floquet  modes,  and  an  integral  over  kx-  Due  to  the  symmetrical  properties  of 
the  Green’s  function,  the  basis  functions  and  the  test  functions,  all  the  integrands 
are  found  to  be  even  functions  of  kx-  So  the  integration  path  can  be  reduced  to 
an  integral  from  0  to  00.  In  numerical  computation,  the  path  of  integration  in  the 
complex  kx  plane  is  deformed  to  avoid  the  singularities  on  the  real  axis  [9]. 

In  the  following  calculations,  the  parameters  used  for  Fig.  2  are:  di  =  62  =  d^  — 
0.2  mm,  p  =  0.5  mm,  wj  =  u>2  =  0.125  mm,  Tj  =  1.7  mm,  and  pj  =  pj,  =  /ip- 
Since  the  crossing  strip  length  is  much  longer  than  the  signal  line  width,  the 
current  near  the  crossing  strip  edges  is  relatively  small  so  that  the  edge  condition 
in  the  basis  functions  can  be  neglected. 
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Figtxre  2.  (a)  Cross  section  of  an  open  structure,  (b)  Cross  section  of  a 

closed  structure. 
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Figure  3  shows  the  dispersion  characteristics  of  a  closed  microstrip  line  loaded 
with  periodic  crossing  strips  (Fig.  2(b)).  The  result  shows  that  the  first  stopband 
appears  due  to  the  coupling  between  the  Floquet  modes  n  =  0  and  n  =  —  1  of 
the  fundsunental  mode  of  the  signal  line.  The  upper  and  lower  bounds  of  the 
stopband  is  denoted  by  otu  and  respectively.  At  higher  frequencies,  higher 
order  stopbands  are  encountered  because  of  the  interaction  with  the  higher  order 
modes  of  the  signal  line.  However,  we  concentrate  only  on  the  first  stopband 
which  is  in  the  region  of  practical  interest. 

In  Fig.  4,  the  dispersion  characteristics  of  an  open  and  a  shielded  structure 
are  plotted  in  solid  and  dashed  lines,  respectively.  It  can  be  seen  that  both  the 
stopband  position  and  width  are  close  to  those  of  each  other.  This  is  because  the 
fields  are  mostly  confined  under  the  first  layer  where  the  coupling  between  signal 
line  and  crossing  strips  takes  place.  So  removing  the  top  conducting  plate  does 
not  affect  the  stopband  properties  much  in  this  case.  This  point  is  illustrated  in 
Figs.  5(a)  and  (b)  which  shows  the  effect  of  changing  di  on  the  stopband  position 
and  width,  respectively.  In  the  following,  we  are  going  to  investigate  the  effect  of 
anisotropy  in  the  second  and  the  third  layer  of  a  closed  structure.  It  is  believed 
that  similar  effects  can  be  observed  in  an  open  structure. 

The  plot  in  Fig.  6  shows  the  effect  of  the  anisotropy  ratio  (AJZ  =  e2/£2z)  of 
the  second  layer  on  the  stopband  position  and  the  stopband  width.  The  center 
frequency  of  the  stopband  is  not  much  affected  by  the  anisotropy.  However,  the 
stopband  width  is  quite  sensitive  to  it.  The  width  increases  with  1/AR.  For  fixed 
<2,  it  corresponds  to  an  increase  of  e2z<  'which  enhances  coupling  between  the 
signal  line  and  the  crossing  strips,  resulting  in  the  rise  of  stopband  width.  For 
1/AR  >  1,  it  is  found  in  the  dispersion  diagram  that  a  high  order  stopband  starts 
to  merge  with  the  first  order  stopband,  resulting  in  a  large  stopband  width.  Fig.  6 
is  thus  plotted  up  to  that  value  only. 

In  Fig.  7,  we  investigate  the  effect  of  anisotropy  in  the  third  layer  on  the 
stopband  properties.  As  we  have  expected,  the  stopband  width  is  not  so  sensitive 
to  the  anisotropy  in  the  third  layer  as  it  is  in  the  second  layer  where  coupling 
occurs.  The  change  of  stopband  position  with  the  anisotropy  is  close  to  that  in 
the  second  layer.  Both  are  due  to  the  change  of  the  dispersion  characteristics  of 
the  signal  line  which  results  in  the  lowering  of  the  intersecting  point  of  the  Floquet 
modes  n  =  0  and  n  =  —1.  A  high  order  stopband  is  encountered  for  1/AR  >  1. 

Various  combinations  of  substrate  materials  have  been  used  to  minimize  the 
stopband  width  for  the  closed  structure  (Fig.  2(b)).  The  results  are  summarized 
in  the  following  table: 


Case 

layer  1 

layer  2 

layer  3 

{kop/ir)c 

A{kop/ir) 

1 

10£0 

sapph. 

lOto 

0.3093 

5.93E-3 

2 

lOto 

lOco 

lOeo 

0.3144 

3.67E-3 

3 

10«o 

Eps-10 

lOfo 

0.3061 

3.40E-3 

4 

2.3«o 

lOeo 

lOeo 

0.3194 

0.80E-3 

5 

2.3eo 

lOfo 

sapph. 

0.3054 

0.53E-3 

6 

«0 

lOeo 

sapph. 

0.3058 

0.87E-3 
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where  (fcoP/’'')c  is  the  normalized  center  frequency  of  the  stopband  and  A(jkoP/w) 
is  the  normalized  stopband  width. 

The  two  types  of  anisotropic  substrates  considered  are  Epsilam-10  (e  =  IScq, 
Cz  =  10.3eo)  &°d  sapphire  (e  =  9.4co,  Cz  =  11.6eo).  Comparison  shows  that  the 
fifth  case  has  the  smallest  stopbsmd  width.  In  fact,  the  stopband  width  is  quite 
sensitive  to  the  separation  of  the  crossing  strips  due  to  the  resonance  effect  [6]. 
Once  the  periodicity  p  is  fixed,  the  stopband  width  can  be  minimized  by  a  proper 
choice  of  substrate  material. 


V.  CONCLUSIONS 

A  dyadic  Green’s  function  formulation  for  the  analysis  of  open  and  closed  noi- 
crostrip  lines  in  the  presence  of  periodic  crossing  strips  in  a  stratified  uniaxially 
anisotropic  medium  is  presented.  The  dispersion  characteristics  for  a  three-layer 
structure  is  studied.  Numerical  results  illustrate  the  relationship  between  the 
stopband  properties  and  the  material  puameters.  The  effect  of  anisotropy  has 
also  been  investigated.  It  is  found  that  the  crossing  strip  separation  and  the 
anisotropy  in  the  second  layer  are  important  factors  affecting  the  stopband  width. 
To  achieve  smaU  stopband  width,  careful  choice  of  anisotropy  must  be  made  to 
avoid  the  lowering  of  the  high  order  stopband.  It  should  also  be  noted  that  by  the 
proper  choice  of  substrate  materials,  the  stopband  width  can  be  much  reduced  for 
fixed  crossing  strip  separation. 
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APPENDIX  A 

Using  the  dyadic  Green’s  function  formulation  [7],  the  elements  of  (ij(fczn)  can 
be  obtained  as  follows: 

For  t  =  j,  where  the  source  and  observation  points  are  in  the  same  layer,  we 
have 

(XI) 

(X2) 

and  for  j  >  i,  where  the  source  is  in  layer  (j)  and  observation  point  is  in  layer  (t), 
we  have 

(X3) 
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where  and  are,  respectively,  the  TE  and  TM  upward  transmission 

coefRcients  from  layer  (j)  to  layer  (i),  given  by 

vTE  _  yTE  ■^u(/+l))  M'S! 

yTM  _  yTM  */  *(I+I)> 

for  I  =  (m  —  2),  (m  —  3), . . . ,  0,  and  for  1  =  (m  —  1),  we  have 
vTE _ (l+-Ru^) 


^U{m— l),m 


fcj  kmt _ 


core 
:  c  I.  L 

L  Ok  C 

c «. 

f  Oi  C 

G  I 

(  O'  : 

-  c  - 


ifcW=  /]fe2_±ijfe2 
ifcW  =  /it?  -  iifci 


*};'  =  y  -  fk]  iAio) 

kj  =  (^11) 

Using  the  symmetrical  properties  of  the  dyadic  Green’s  function  in  the  layered 
media,  it  can  be  shown  that  /jj  =  /,“,  where  a  denotes  TE  or  TM.  In  the  above 
equations,  the  superscripts  (h)  and  (e)  denote  TE  and  TM  fields,  respectively. 
R^l  and  iZJ/  are,  respectively,  the  Fresnel  reflection  coefficients  at  the  lower  and 
upper  boundaries  of  layer  (I)  and  can  be  determined  recursively  by  the  following 
relations 


1  +  i?“  R?  e*'*r'+i)*‘*'+i 


(A12) 


where  I  =  0, . . .  ,(n  —  1)  and  J2^„  =  IJJj. 
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(A13) 


where  /  =  2, 3, . . . ,  i  and  R^j  =  •Rj'o-  ■^j_i)  '^J(J+1)  Fresnel  reflection 

coeffldents  across  the  interface  between  layers  (?)  and  (1  +  1). 
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Abstract 


Problem  Configuration  and  Model 


The  electromagnetic  radiation  from  a  VLSI  chip  package 
and  heatsink  structure  is  analysed  by  means  of  the  finite-dif¬ 
ference  time-domain  (FD-TD)  method.  The  FD-TD  algorithm 
implemented  incorporates  a  multi-zone  gridding  scheme  to  ac¬ 
commodate  fine  grid  cells  in  the  vicinity  of  the  heatsink  and 
package  cavity  and  sparse  gridding  in  the  remainder  of  the  com¬ 
putational  domain.  The  issues  pertaining  to  the  effects  of  the 
heatsink  in  influencing  the  overall  radiating  capacity  of  the  con¬ 
figuration  are  addressed.  Analyses  are  facilitated  by  using  sim¬ 
plified  heatsink  models  and  by  using  dipole  elements  as  sources 
of  electromagnetic  energy  to  model  the  VLSI  chip.  The  po¬ 
tential  for  enhancement  of  spurious  emissions  by  the  heatsink 
structure  is  illustrated.  For  heatsinks  of  typical  dimensions, 
resonance  is  possible  within  the  low  gigahertz  frequency  range. 
The  potential  exploitation  of  the  heatsink  as  an  emissions  shield 
by  appropriate  implementation  schemes  is  discussed  and  evalu¬ 
ated. 


Introduction 

The  extension  of  regulatory  electromagnetic  emissions  test 
limits  to  the  gigahertz  range  of  frequencies  for  high  performance 
computing  equipment  has  necessitated  caution  in  the  design  and 
implementation  of  components  down  to  the  VLSI  chip  package 
level  for  computing  systems.  At  issue  is  the  avoidance  of  poten¬ 
tially  efficient  radiators  at  the  design  and  development  phase. 
The  objective  is  to  minimize  the  need  for  cost  prohibitive  cor¬ 
rective  measures  commonly  invoked  in  reaction  to  problems  en¬ 
countered  in  tests  on  assembled  systems. 

In  this  paper,  the  effects  of  a  heatsink  over  a  chip  package 
on  the  electromagnetic  emissions  characteristics  of  the  chip  and 
package  environment  are  addressed.  Such  an  interest  is  justified 
by  the  following  reasons.  First,  the  heatsink  is  typicaUy  metal¬ 
lic  and  is  in  close  proximity  to  a  major  energy  source  —  the 
chip.  Also,  the  electrical  dimensions  of  the  structure  are  com¬ 
parable  to  the  wavelengths  of  concern.  Finally,  in  light  of  the 
ever-increasing  power  levels  of  high  performance  integrated  cir¬ 
cuitry,  the  heatsink  is  an  indispensable  component  of  the  VLSI 
chip  packaging  configuration.  The  objectives  are  to  identify  and 
quantify  the  effects  of  the  heatsink  on  the  radiating  properties  of 
the  package  structure;  to  understand  and  exploit  the  radiation 
mechanisms;  and  to  evaluate  viable  heatsink  implementation 
schemes  for  minimizing  the  overall  radiating  capacity. 


•  Current  .ddr««:  WAVETRACER,  Acton.  MA  01720. 


The  electromagnetic  radiation  properties  of  a  heatsink  and 
integrated  circuit  (IC)  package  configuration  are  analysed  by 
means  of  suitable  models.  A  typical  heatsink/package  config¬ 
uration  of  interest  is  shown  in  Figure  1.  The  heatsink  may  be 
modeled  as  a  perfectly  conducting  rectangular  slab  positioned 
over  a  finite-size  dielectric  medium  representing  the  chip  pack¬ 
age.  The  chip  is  supported  in  turn  by  a  dielectric  layer  of  infinite 
extent  over  an  infinite  ground  plane  which  models  the  substrate 
or  printed  wiring  board  (PWB)  with  (at  least)  one  reference 
layer.  The  electromagnetic  sources  used  to  model  the  active 
chip  include  electric  and  magnetic  dipoles  of  verticsJ  and  hor¬ 
izontal  orientation  positioned  on  the  package/substrate  dielec¬ 
tric  interface.  The  resulting  model  of  an  electromagnetically- 
coupled  slab  structure  is  shown  in  Figure  2.  The  generation  of 
a  model  with  reduced  complexity  for  both  source  and  heatsink 
is  to  facilitate  the  interpretation  of  results  based  on  underlying 
physics.  The  extension  to  more  complex  heatsink  (e.g.,  finned) 
and  source  models  may  be  readily  accomplished  with  the  rect¬ 
angular  gridding  scheme  associated  with  the  numerical  solution 
employed.  With  sufficiently  fine  grids,  circular  heatsinks  may 
also  be  adequately  modeled.  In  addition,  the  models  are  modi¬ 
fied  to  to  reflect  alternate  implementation  schemes  for  emissions 
suppression  purposes  in  order  to  conduct  quantitative  analyses. 
Field  strengths  and  radiated  powers  are  computed  for  quanti¬ 
tative  analysis  and  evaluation. 


Figure  1.  Cutaway  view  of  a  VLSI  package/heatsink  configu¬ 
ration. 


CH3044-5/9 1/0000-0393  S01.00  ®  1991  IEEE 


393 


t 


9«oufid  pltftt 


Figure  2.  Cross-sectional  view  of  simplified  package/heatsink 
model. 

Method  of  Solution 

The  method  employed  in  the  analysis  is  the  finite-difference 
time-domain  (FD-TD)  technique.  It  is  based  on  the  discretiza¬ 
tion  of  the  electric  and  magnetic  fields  over  rectangular  grids 
together  with  the  finite  difference  approximation  of  the  spa¬ 
tial  and  temporal  derivatives  appearing  in  the  differential  form 
of  Maxwell's  equations.  The  reasons  for  which  the  FD-TD 
methodology  was  selected  include  the  relative  ease  of  imple¬ 
mentation  for  complicated  geometries,  the  requirement  of  only 
simple  arithmetic  operations  in  the  solution  process,  and  the 
flexibility  for  time-  and  frequency-domain  analyses. 


electric  and  magnetic  fields.  To  achieve  accurate  results,  the 
cell  sizes  are  taken  to  be  a  fraction  of  the  smallest  wavelength. 
The  time  increment  and  the  cell  size  are  related  by  the  stability 
criterion  (2], 


where  c,  is  the  speed  of  light  in  free  space.  Fields  are  set  to  be 
zero  initially  everywhere  to  satisiy  the  causality  condition  con¬ 
sistent  with  zero  excitation  for  time  less  than  zero.  The  bound¬ 
ary  conditions  are  continuity  of  tangential  electric  and  magnetic 
fields  on  material  interfaces,  vanishing  tangential  electric  fields 
on  perfect  conductors,  and  the  absorbing  boundary  conditions 
on  the  boundary  of  the  computational  domain.  The  absorb¬ 
ing  boundary  conditions  [3]  are  used  to  limit  the  computational 
domain  by  simulating  unbounded  space.  The  minimum  dis¬ 
tance  to  the  absorbing  boundary,  i.e.,  computational  boundary, 
from  the  heatsink  is  deternained  by  consideration  of  reflection 
error,  computation  time,  and  memory.  For  a  given  cell  size,  the 
reflection  error  decreases  but  the  computation  time  increases 
substantially  with  increasing  distance  to  the  boundary. 


/'(l/Az)’-l-(l/Ay)*-f-(l/Az)> 


In  the  FD-TD  technique,  a  computational  domain  is  first 
defined  and  divided  into  rectangular  cells.  Electric  and  mag¬ 
netic  fields  are  spatially  discretized  in  a  staggered  manner  [1] 
SIS  shown  in'  Figure  3.  Electric  fields  are  assigned  to  half-integer 
(n  -H  1/2)  time  steps  and  magnetic  fields  are  assigned  to  integer 
(n)  time  steps  for  the  temporal  discretization  of  fields.  Next, 
the  spatial  and  temporal  derivatives  of  the  two  Maxwell's  curl 
equations  are  approximated  using  center  differences.  Maxwell's 
curl  equations  for  a  time-  and  frequency-invariant  medium  are: 

«  77  dE  „  BE 

V  y.  H  =  a,E -i- V  x  £  = 

where  <o  is  the  free-space  permittivity,  8.854  X  10“**  F/m  and 
fie  is  the  free-space  permeability,  4ir  x  10“*  H/m.  In  addition, 
t,  and  fir  are  respectively  the  relative  permittivity  (dielectric 
constant)  and  relative  permeability  of  the  medium;  while  ff, 
is  the  electric  conductivity.  Maxwell's  divergence  equations  are 
ignored  since  the  curl  equations  with  appropriate  boundary  con¬ 
ditions  uniquely  determine  the  solution.  In  rectilinear  coordi¬ 
nates,  the  curl  equations  are  rewritten  as: 

BH.  BH.  r-  BEy  BH„ 

By  8z  Bt  ay  Oz  at 

BH^  BH.  r-  9Ey  BE.  BE.  BH. 

Bz  ~  Bz  +  gi  '  g^  = -/'r/i, 

BH.  BH.  „  BE.  BE.  BE.  BH. 

Bz  By  gt  <  g^  g^  = ^ 

Difference  equations  are  derived  from  these  six  equations  by  ap¬ 
plying  center  differencing.  For  example,  the  expression  relating 
H.  to  E.  and  E,  at  time  (n  -f  1/2) At  is  given  by: 


Ay 


The  center  difference  ensures  that  the  spatial  and  temporal  dis¬ 
cretizations  are  o!  second  order  where  errors  are  proportional 
to  the  square  of  the  cell  size  and  time  increment.  Finally,  with 
appropriate  initial  and  boundary  conditions,  the  solutions  to 
the  difference  equations  are  obtained  through  explicit  leapfrog 
time  marching.  This  corresponds  to  alternating  the  adv»nce  of 
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Figure  3.  Staggered  rectangular  grid  on  typical  unit  cell  in 
computational  domain. 

The  size  of  the  discretization  cells  should  be  reasonably 
small  in  order  to  model  features  of  the  heatsink  and  capture 
correct  field  variations.  On  the  other  hand,  since  the  heatsink 
should  not  be  too  close  to  the  computational  boundary,  the  cell 
count  is  highest  outside  the  heatsink  when  a  small  cell  size  is 
used.  To  avoid  the  excessive  number  of  cells  which  would  other¬ 
wise  be  generated,  a  three-zone  gridding  scheme  [4j  is  used  (Fig¬ 
ure  4).  The  first  zone  contains  the  heatsink  and  dipole  source; 
it  has  the  finest  grid  among  the  three  zones.  The  second  zone 
has  the  same  height  as  the  first  and  extends  horizontally  from 
the  heatsink  to  the  outer  boundary  of  the  computational  do¬ 
main.  The  cell  size  of  this  zone  has  the  same  vertical  dimension 
as  those  of  the  first  zone;  while  the  horizontal  dimensions  may 
be  a  few  times  larger  than  the  first  zone.  The  remainder  of  the 
computational  domain  belongs  to  the  third  zone.  The  horizon¬ 
tal  dimensions  of  the  cells  are  identical  to  those  of  the  second 
zone,  and  the  vertical  dimension  of  the  cells  may  be  a  few  times 
larger  than  the  first  or  second  zones.  Fields  at  nodes  on  an  in¬ 
terface  between  two  zones  are  calculated  via  a  combination  of 
parabolic  curve  fitting  and  linear  interpolation.  For  example,  in 
Figure  5,  the  partial  derivative  of  Hy  with  respect  to  the  ver¬ 
tical  direction  z  at  position  of  at  the  nth  time  step  if 

given  by 
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Figure  4.  A  three-rone  gridding  scheme  for  multi-rone  grid- 
ding.  The  heatsink  is  in  rone  1. 

8  ffn 

dz  ^  A{\  +  a){3  +  a)  »<■•>•*> 
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Second  order  accuracy  is  therefore  maintained  everywhere.  Al¬ 
though  the  discretisation  scheme  is  second  order  everywhere, 
the  proportionality  constants  for  the  discretisation  errors  are 
dilTerent  due  to  different  cell  sires.  Reasonable  accuracy  it  main¬ 
tained  due  to  the  fact  that  the  spatial  variation  of  the  helds 
outside  the  first  rone  is  relatively  small  and  that  all  cell  sizes 
remain  a  fraction  of  the  smallest  wavelength.  The  size  of  the 
time  step  is  primarily  determined  by  the  smallest  grid  size.  For 
simplicity,  the  variable  time  step  implementation  [5]  is  not  used. 

The  fields  beyond  the  computational  domain  are  calculated 
using  Huygens'  principle  [6].  The  post-processing  of  full  time 
domain  information  for  the  purpose  of  generating  frequency  do¬ 
main  data  is  expensive  in  storage  and  is  inefficient.  The  fields 
are  instead  calculated  at  selected  multiple  frequencies.  To  ob¬ 
tain  multi-frequency  data,  the  time  waveforms  of  the  dipole 
excitations  are  chosen  to  be  modulated  Gaussian  pulses,  i.e., 
~  cos(w,(t  -  to))  expl-((t  -  t,)/r)*l,  where  w.  is  the  center  fre¬ 
quency  of  the  excitation,  t,  is  the  delay,  and  T  determines  the 
pulse-width.  Complex  amplitudes  at  each  frequency  are  calcu¬ 
lated  simultaneously  on  a  selected  surface  outside  the  first  zone 
using  the  discrete  Fourier  transform. 


Results  and  Discussion 

The  basic  model  described  in  Figure  2  is  analysed  with 
the  use  of  horizontal  electric  and  magnetic  dipoles  and  vertical 
electric  and  magnetic  dipoles  as  sources.  The  choice  of  such 
elements  circumvents  the  problem  of  deriving  a  rigorous  model 
for  the  VLSI  circuitry  as  an  electromagnetic  source  and  elim¬ 
inates  factors  which  may  otherwise  clutter  the  understanding 
and  interpretation  of  the  physical  processes  associated  with  the 
spurious  radiation  study.  Of  these  sources  only  the  horizon¬ 
tal  magnetic  dipole  (HMD)  and  vertical  electric  dipole  (VED) 
show  radiation  enhancement  in  the  presence  of  the  heatsink. 
The  frequency  range  considered  does  not  exceed  10  GHz.  The 
dominant  horizontal  (electric  field)  polarization  generated  by 
the  other  2  sources  remains  below  cutoff  within  this  frequency 


Figure  5.  Node  and  field  designations  for  multi-zone  gridding 
in  vertical  direction. 

range  and  for  the  typical  dimensions  considered.  Both  classes 
of  sources  may  be  found  in  VLSI  package  configurations  —  the 
VED  in  the  form  of  vias  and  pins  and  the  HMD  as  vertically 
stratified  signal  and  return  paths.  Results  and  discussions  will 
focus  on  one  or  both  of  the  VED,  HMD  sources.  Unless  stated 
otherwise,  the  dimensions  of  the  conducting  slab  are  4.8  cm 
square  and  2.5  mm  thick.  Its  lower  surface  is  located  7.5  mm 
above  a  lower  ground  or  referenee  plane.  In  the  results  pre¬ 
sented,  the  effects  of  the  dielectric  layers  are  not  considered. 
The  dipole  sources  are  located  at  the  nodes  closest  to  the  cen¬ 
ter  of  the  heatsink-ground  plane  cavity. 


Figure  6  compares  the  normalized  radiated  power  of  an 
HMD  over  a  ground  plane  in  the  presence  and  absence  of  a 
heatsink  for  a  frequency  range  of  1  to  6.5  GHz.  Normalization 
is  with  respect  to  the  total  radiated  power  of  an  identical  dipole 
in  unbounded  space.  The  resonant  behavior  around  2.4  and  5-5 
GHz  in  the  presence  of  the  heatsink  is  evident,  with  the  first  res¬ 
onance  showing  a  significant  enhancement  of  radiation  over  the 
case  without  the  heatsink.  This  resonance  feature  is  illustrated 
in  Figure  7  where  the  time  response  of  the  vertical  Ei-field  at 
the  indicated  location  on  the  periphery  of  the  heatsink-ground 
plane  cavity  is  shown.  The  oscillations  at  the  resonant  frequency 
become  increasingly  evident  after  700  time-steps.  The  quality 
factor  (Q)  of  the  structure  may  also  be  deduced  from  the  decay 
rate  and  the  resonant  frequency.  Figure  8  shows  the  surface 
plots  of  the  vertical  electric  field  amplitudes  over  a  horizontal 
plane  that  traverses  the  heatsink  and  ground  plane  cavity  and 
spanning  the  computational  domain  at  4  different  time  steps. 
Each  plot  is  normalized  to  the  maximum  value  for  the  respective 
time  step.  The  absence  of  any  obvious  reflected  disturbances 
from  the  edges  of  the  computational  domain  illustrates  the  ef¬ 
ficiency  of  the  absorbing  boundary  condition.  The  sustained 
resonant  mode  pattern  is  evident  after  1000  time  steps. 

Figure  9  shows  the  corresponding  results  over  the  frequency 
range  of  1  to  8  GHz  with  a  VED  as  the  source.  Again,  substan¬ 
tial  enhancement  in  radiated  power  in  the  presence  of  a  heatsink 
occurs  at  the  primary  resonance  of  4.8  GHz,  and  at  the  second 
resonance  near  7  GHz.  The  symmetry  and  anti-symmetry  of  the 
dipole  source  fields  are  responsible  for  the  difference  in  excited 
resonances  with  the  HMD  and  VED  dipoles. 
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FRfQUENCY  (CHi) 

Figure  6.  Comparison  of  normalized  total  radiated  power  in 
presence  and  absence  of  a  heatsink  with  a  HMD  source.  Source 
function  ~  cos(u,(t-to))exp[-((t-l»)/r)’l,  where  u»o  =  Ssr  x 
10*  rad/s,  =  4T  =  0.8  ns. 
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TIME  STEPS 

Figure  7.  Time  response  of  vertical  electric  field  at  indicated 
position  on  the  boundary  of  heatsink  and  ground  plane  cav¬ 
ity  for  HMD  source  of  Figure  6.  Each  time  step  corresponds 
approximately  to  2.2  picoseconds,  fn  s  2.3  GHz. 
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n  s  1000 


Figure  8.  Vertical  electric  field  amplitude  plots  over  horizontal 
plane  through  heatsink/ground  plane  cavity  for  HMD  source  of 
Figure  6.  n  is  the  number  of  time  steps  elapsed,  and  each  time 
step  corresponds  approiumately  to  2.2  picoseconds. 
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Figure  0.  Comparison  of  normalized  total  radiated  power  in 
presence  and  absence  of  heatsink  with  an  VED  source.  Source 
function  ~  cos(a)e(t  —  ta))ejcp(— ((<  —  <,)/r)’),  where  w#  =  5irx 
10*  rad/s,  t,  =  AT  =  0.8  ns. 


Although  modeled  so  far  as  a  suspended  slab  of  conduct¬ 
ing  material,  the  heatsink  may  also  be  supported  by  posts  (e.g., 
Figure  1)  which  may  or  may  not  be  connected  to  any  other  con¬ 
ducting  paths.  In  Figure  10  we  show  the  model  of  a  heatsink 
supported  with  4  posts,  of  4  mm  square  cross-section,  connected 
to  a  reference  or  ground  plane.  We  also  investigate  the  effect  of 
increasing  the  number  of  the  posts  as  a  means  of  exploiting  the 
heatsink  as  an  emissions  shield  through  appropriate  implemen¬ 
tation.  This  is  particularly  pertinent  in  light  of  the  observation 
of  radiation  enhancement  by  the  indispensable  heatsink.  The 
rationale  for  consideration  of  such  posts  is  the  reduction  in  the 
size  of  the  wavelengths  allowed  to  leak  through  the  smaller  gaps, 
whose  effectiveness  is  weakened  by  the  interaction  between  the 
resulting  multiple  gaps. 
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Figure  10.  Model  for  heatsink  with  4  support  posts  connected 
to  reference  plane. 


Figure  11  shows  the  effects  of  the  implementation  of  4  a; 
8  uniformly  distributed  posts  on  the  radiated  power.  The  dipt 
source  is  a  HMD.  The  figure  shows  the  typical  effects  of  addi: 
posts  —  the  shifting  of  resonant  frequencies  upwards  which  al 
renders  the  resulting  configuration  to  be  of  lower  Q.  In  impi 
menting  such  an  option  for  emissions  suppression,  care  has  to  ) 
taken  to  ensure  that  such  shifts  move  the  resonant  frequenci 
beyond  the  range  of  test  frequencies. 


Figure  11.  Comparison  of  normalized  total  radiated  powe 
with  HMD  source  for  heatsink  with  4  and  8  posts  and  withou 
posts. 


Another  potential  means  of  exploiting  the  heatsink  foi 
emissions  suppression  is  to  employ  a  skirt  of  conducting  mate 
rial  beneath  and  around  the  periphery  of  the  heatsink,  therebj 
serving  as  a  compliant  containment  which  would  dissipate  thi 
contained  energy  with  its  finite  conductivity.  Figure  12  shows 
the  model  considered.  The  gaps  in  the  gasket  material  are  tc 
facilitate  lead  passage.  The  simulation  results,  comparing  radi¬ 
ation  from  a  gasketed  and  ungasketed  heatsink  with  an  HMD 
as  source,  are  shown  in  Figure  13  where  a  conductivity  of  10 
mhos/m,  a  dielectric  constant  of  4,  and  a  thickness  of  approxi¬ 
mately  7.5  mm  are  assumed  for  the  gasket  material.  This  con¬ 
ductivity  is  conservative  since  typical  bulk  conductivities  for 
conductive  gasket  material  are  of  the  order  of  10*  mhos/m.  De¬ 
spite  the  conservative  estimate  on  the  conductivity  of  a  typi¬ 
cal  material,  the  improvement  in  emissions  reduction  near  reso¬ 
nance  it  substantial,  although  the  lower  Q  that  results  may  tend 
to  raise  emission  levels  at  nearby  frequencies.  For  reference,  the 
size  of  the  skin  depth  at  1  GHz  for  the  conductivity  assumed  is 
approximately  5  mm  and  will  decrease  as  the  square-root  of  the 
frequency.  The  frequency  response  indicaf-s  a  sb'ght  downward 
shift  in  the  resonant  frequency  due  to  an  •*  fectively  denser  di¬ 
electric  medium  beneath  the  heatsink  and  a  broader  bandwidth 
corresponding  to  a  lowered  Q  due  to  added  energy  dissipation 
in  the  gasket  material. 
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Figure  12.  Model  for  heatsink  with  a  conducting  gasket. 


FRfgUENCY  (CUt) 


Figure  13.  Comparison  of  normalized  total  radiated  power 
with  HMD  source  for  heatsink  with  and  without  gasketing. 


For  the  purpose  of  suppressing  emissions,  exclusive  usage  of 
posts  appears  viable  if  the  resonant  frequencies  may  be  shifted 
beyond  the  frequency  range  of  concern.  The  use  of  gaskets 
results  in  slightly  lower  resonant  frequencies  with  broader  fre¬ 
quency  responses.  This  spreads  the  resonance  effect  over  a 
broader  frequency  range,  albeit  at  lower  amplitudes.  A  rec- 
onunendation  would  be  to  use  a  combination  of  the  gasketing 
and  multiple  posts  schemes.  The  shielding  performance  of  con¬ 
ductive  gaskets  are  expected  to  improve  at  higher  frequencies  to 
which  multiple  posts  shift  the  resonant  frequencies  of  concern. 


The  FD-TD  techm'que  has  been  applied  to  the  analyses  of 
the  radiation  from  a  VLSI  heatsink  and  package  configuration. 
The  technique  was  chosen  for  its  flexibility  in  treating  this  class 
of  problems  involving  potentially  complex  configurations  where 
neither  approximate  nor  analytical  methods  are  practical.  The 
multi-zone  gridding  scheme  implemented  allowed  a  high  grid 
resolution  in  the  vicinity  of  the  heatsink  without  sacrificing  over¬ 
all  cell  number  thereby  allowing  increased  numerical  efficiency 
for  low  frequencies.  Additional  advantages  of  this  discretization 
scheme  include  improved  geometrical  and  material  modeling. 
The  analysis  of  the  heatsink,  a  critical  component  within  the 
computer  packaging  environment,  illustrated  the  significance  of 
resonance  due  to  appreciable  electrical  dimensions  on  the  spu¬ 
rious  electromagnetic  radiation  from  the  VLSI  packaging  envi¬ 
ronment.  Simulations  with  typical  heatsinks  dimensions  showed 
the  occurrence  of  resonance  in  the  low  gigahertz  range.  The  ef¬ 
fects  of  the  presence  of  a  heatsink  on  the  radiation  properties 
of  dipole  models  have  been  explained  and  the  features  of  res¬ 
onant  behavior  presented.  The  effectiveness  and  consequences 
of  exploiting  practical  heatsink  implementation  options  such  as 
grounding  and  shielding  to  reduce  electromagnetic  emissions  has 
been  discussed  and  quantified  using  the  FD-TD  numerical  tech¬ 
nique. 
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Because  the  effects  of  diffraction  during  proximity-print  x-ray 
lithography  are  of  critical  importance,  a  number  of  previous 
researchers  have  attempted  to  calculate  the  diffraction  patterns  and 
minimum  achievable  feature  sizes  as  a  function  of  wavelength  and 
gap.  Work  to  date  has  assumed  that  scalar  diffraction  theory  is 
applicable — as  calculated,  for  example,  by  the  Rayleigh-Sommerfeld 
formulation — and  that  Kirchhoff  boundary  conditions  can  be  applied. 
Kirchhoff  boundary  conditions  assume  that  the  fields  (amplitude  and 
phase)  are  constant  in  the  open  regions  between  absorbers,  and  a 
different  constant  in  regions  just  under  the  absorbers  (i.e.,  that  there  are 
no  fringing  fields).  An  x-ray  absorber  is,  however,  best  described  as  a 
lossy  dielectric  that  is  tens  or  hundreds  of  wavelengths  tall,  and  hence 
Kirchhoff  boundary  conditions  are  unsuitable.  In  this  report  we  use 
two  numerical  techniques  to  calculate  (on  a  Cray  2  supercomputer) 
accurate  diffracted  fields  from  gold  absorbers  for  two  cases:  a  30  nm- 
wide  line  at  X  =  4.5  nm,  and  a  100  nm-wide  line  at  X  =  1.3  nm.  We 
show  that  the  use  of  Kirchhoff  boundary  conditions  introduces 
unphysically  high  spatial  frequencies  into  the  diffracted  fields.  The 
suppression  of  these  frequencies — which  occurs  naturally  without  the 
need  to  introduce  an  extended  source  or  broad  spectrum — improves 
exposure  latitude  for  mask  features  near  0.1  pm  and  below. 
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I.  INTRODUCTION 


The  limits  and  practicality  of  proximity-print  x-ray  lithography  continues  to 
be  a  topic  of  discussion  and  debate.  Of  particular  concern  are  the  linuts  of 
resolution  imposed  by  the  effects  of  diffraction.  Because  the  mechanical 
limits  imposed  on  the  mask-substrate  gap  during  the  volume  manufacturing 
of  ULSI  circuits  are  not  presently  certain,  it  remains  prudent  to  ask:  what  is 
the  minimum  practical  feature  size  that  can  be  printed  at  a  given  gap? 

In  order  to  resolve  these  issues,  a  series  of  papers  have  appeared  in  the 
literature  which  present  theoretical  calculations  for  the  diffraction  of  x  rays 
from  mask  absorbers.  Early  papers  simply  considered  absorption  [1-3],  but 
later  papers  also  included  the  effects  of  phase-shift  [4-8].  Most  recently, 
authors  have  included  the  effects  of  source  spatial  and  temporal  coherence, 
and  have  generated  exposure  "trees"  which  allow  the  determination  of 
exposure  latitude  versus  gap  for  various  types  of  mask  features  [9-15]. 

The  method  most  commonly  used  to  calculate  the  diffraction  pattern  is  to 
apply  a  Fresnel-Kirchhoff  or  Rayleigh-Sommerfeld  diffraction  integral  [16,17], 
or  a  more  sophisticated  formulation  (Hopkin's  formula)  which  takes  into 
accoimt  source  partial  coherence  [18].  These  calculations  can  be  carried  out 
either  in  the  spatial  or  the  spatial-frequency  domain.  In  any  of  these  cases, 
approximate  boundary  conditions  known  as  Kirchhoff  boundary  conditions 
(KBC)  are  generally  applied.  KBC  assume  that  the  field  (amplitude  and  phase) 
is  constant  in  the  open  regions  between  absorbers,  and  also  constant  (but 
attenuated  and  phase  shifted)  in  regions  just  under  the  absorbers — in  other 
words,  that  there  are  no  fringing  fields  at  the  boundary  between  the  two 
regions. 

In  general,  KBC  apply  when  the  wavelength,  X,  is  much  smaller  than  the 
lateral  size,  d,  of  the  feature  being  printed — which  is  the  case  in  most  x-ray 
lithography  (e.g.,  d  =  0.1  fim,  X  =  1  nm).  However,  what  is  not  generally 
recognized  is  that  KBC  will  not  necessarily  apply  when  the  absorbers  are  lossy 
dielectrics  that  are  tens  or  hundreds  of  wavelengths  tall. 

Work  to  date  concerning  the  printability  of  0.25  4m  features  is  probably 
reasonably  accurate.  However,  in  this  report  we  show  that  the  assumption  of 
KBC  for  features  near  0.1  pm  and  below  is  not  tenable.  In  particular,  we  show 
that  the  suppression  of  the  undesirable  high-spatial-frequency  components — 


2 


which  some  authors  note  an  extended  source  and/or  a  broad  spectrum 
achieves — occurs  nahirally  in  the  absorber  due  to  the  "lossy  dielectric"  effect. 

IL  CALCULATIONS 

We  used  two  different  methods  to  calculate  the  diffraction  from  gold 
absorbers:  the  Method  of  Moments  (MoM)  and  the  Finite  Difference-Time 
Domain  (FD-TD)  method.  (We  used  two  algorithms  in  order  to  check  for 
consistency.)  The  only  approximations  inherent  in  these  methods  are  in  the 
discretization  of  the  object  space  and  Maxwell's  equations.  Because  the 
discretization  can  performed  on  a  scale  that  is  small  compared  to  the 
wavelength,  and  furthermore  the  discretization  scale  can  be  reduced  until 
convergence  is  achieved,  accuracy  is  assured. 

Of  these  two  techniques,  the  MoM  is  typically  faster  and  uses  less  memory 
for  single-frequency  calculations.  On  the  other  hand,  the  FD-TD  method  is 
simpler  to  code  and  therefore  less  likely  to  have  errors.  In  practice  we  ran 
small  test  cases  using  both  methods  and  then  increased  the  spatial  resolution 
(reduced  discretization  scale)  until  both  methods  converged  to  the  same 
solution.  Then  the  computationally-intensive  cases  reported  here  were 
calculated  with  the  MoM. 

A.  METHOD  OF  MOMENTS 

The  Method  of  Moments  (MoM)  is  a  numerical  technique  useful  in  the 
solution  of  steady-state  electromagnetic  wave  scattering  and  radiation 
problems  [19,20].  The  method  calculates  steady-state  fields  on  the  surface  of  a 
closed  dielectric  object  in  free  space,  which  is  impinged  upon  by  a  known 
exciting  wave.  The  surface  of  the  object  is  broken  up  into  small  patches 
which  are  small  compared  to  the  wavelength.  The  surface  currents  at  each 
patch  and  thus  the  tangenticil  surface  fields  are  then  calctilated.  Computation 
time  was  up  to  two  hours  on  a  Cray  2.  Once  the  fields  are  known  on  the 
boimdaiy  of  the  object  they  can  readily  be  propagated  to  any  desired  point  or 
plane  [16,17]. 
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B.  FINITE  DIFFERENCE-TIME  DOMAIN  METHOD 


The  Finite  Difference-Time  Domain  (FD-TD)  method  is  a  numerical 
technique  useful  in  the  solution  of  time-dependent  electromagnetic  wave 
scattering  and  radiation  problems  [21-24].  The  method  involves  the 
formation  of  a  computational  domain  which  encompasses  the  object  of 
interest  and  is  typically  several  times  larger.  The  entire  domain  inside  the 
boundary — including  the  object — is  discretized  on  a  rectangular  grid.  The 
spacing  between  adjacent  nodes  on  the  grid  is  small  compared  to  the 
wavelength.  A  clock  is  started  and  incremented  in  time  steps  that  are  small 
compared  to  the  light-travel  time  between  adjacent  nodes.  Then  discretized 
forms  of  Maxwell's  equations  are  used  to  calculate  the  fields  at  each  node 
from  the  fields  at  nearby  nodes  which  were  in  effect  at  the  previous  time  step. 
Absorbing  boundary  conditions  are  imposed  at  the  edges  of  the 
computational  domain  in  order  to  simulate  unbounded  space.  Also,  a 
boundary  condition  is  typically  imposed  on  a  surface  surrounding  the  object 
to  simulate  an  incoming  plane  wave.  The  result  can  be  displayed  as  a  video 
image  of  the  fields  inside  the  domain.  Computation  time  was  up  to  one  hour 
on  a  Cray  2.  In  practice  the  calculation  is  run  until  steady  state  is  achieved, 
and  the  fields  at  nodes  along  a  line  just  under  the  absorber  are  saved  for 
comparison  with  the  MoM  results. 

III.  RESULTS 

Calculations  were  performed  using  MoM  and  FD-TD  on  single  gold 
parallelepiped  absorbers,  infinite  in  length  and  rectangular  in  cross-section, 
which  were  impinged  upon  by  a  monochromatic  plane  wave.  The  electric- 
field  polarization  was  used  (E-field  perpendicular  to  the  page).  (The 
magnetic-field  polarization  yielded  results  similar  to  the  electric-field 
polarization.)  These  were  compared  to  the  results  of  a  Rayleigh-Spmmerfeld- 
Kirchhoff  (RSK)  calculation  [16].  We  considered  two  cases  (see  Table  I),  which 
were  selected  in  order  to  explore  a  range  of  spatial  frequencies:  Case  1,  which 
is  a  30x100  nm  absorber  with  the  4.5  nm  (Cr)  x  ray,  and  Case  2,  which  is  a 
100x250  nm  absorber  with  the  1.3  nm  (Cul)  x  ray.  Note  that  the  attenuation  of 
the  absorber  in  both  cases  is  roughly  equivalent.  (~12  dB). 
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A.  CASE  1  a  =  4.5  run) 

Hgs.  1  (a)  to  (d)  show  the  resulting  intensities  or  irradiance  distributions 
for  both  the  MoM  and  RSK  methods  calculated  for  a  from  0-1.5  (gap  G  from  0- 
0.3  Jim).  Here  a  is  a  dimensionless  gap  given  by  a  =  GX/W^,  where  G  is  the 
gap,  X  is  the  wavelength,  and  W  is  the  linewidth.  Note  the  suppression  of 
high  spatial  frequencies  in  the  MoM  calculation,  and  the  "fuzzy  edge"  of  a  few 
tens  of  nanometers  in  extent.  Figs.  2  (a)  and  (b)  show  the  beneficial  effects  of 
the  suppression  of  high  spatial  frequencies  on  exposure  latitude  in  the  form 
of  exposure  "trees"  [9-12]  which  plot  ±10  %  linewidth  contours  versus  a.  Here 
we  have  used  a  line  bias  of  33  %  (40  nm  resist  line).  (A  line  bias  is  the  use  of  a 
smaller-than-desired  feature  size  in  the  mask  than  on  the  wafer  in  order  to 
compensate  for  diffractive  broadening.)  An  enlightening  way  to  view  the 
suppression  of  high  spatial  frequencies  is  to  plot  the  intensity  of  the  waves  at 
a  =  0  (gap  G  =  0)  as  a  fimction  of  spatial  frequency,  as  shown  in  Fig.  3.  Note 
that  at  low  spatial  frequencies  the  RSK  and  MoM  calculations  agree,  but  the 
MoM  calculation  "rolls  off"  at  around  0.05  nm*l. 

B.  CASE  2  (X  =  1.3nm) 

The  results  for  the  intensities  versus  a  for  Case  2  are  shown  in  Figs.  4  (a)  to 
(d).  Note  that  even  though  the  wavelength  is  much  smaller,  the  "fuzzy  edge" 
length  is  roughly  the  same — a  few  tens  of  nanometers.  This  may  be  due  to 
the  smaller  wavelength  (3.3x  smaller)  being  partially  compensated  by  a  taller 
absorber  (2.5x  taller).  The  exposure  trees  are  shown  in  Figs.  5  (a)  and  (b).  Here 
we  have  used  a  line  bias  of  50%  (150  nm  resist  line).  The  intensity  versus 
spatial  frequency  is  shown  in  Fig.  6.  Note  that  the  roll-off  in  this  case  is  still 
around  0.05  nm'^,  but  that  this  represents  a  higher  spatial  frequency  relative 
to  the  information  content  in  the  larger-line/smaller-wavelength  case. 

IV.  CONCLUSION 

We  have  shown  that  the  use  of  Kirchhoff  boundary  conditions  introduces 
unphysically  high  spatial  frequencies  into  the  diffracted  fields.  The  natural 
suppression  of  these  frequencies  by  the  electromagnetic  properties  of  x-ray 
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absorbers  tremendously  improves  exposure  latitude  for  mask  features  near 
0.1  pm  and  below. 
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Figure  Captions 

FIG.  1-  Intensity  vs.  lateral  position  for  Case  1.  (a)  a  =  0  (gap  =  0),  (b)  a  =  0.5 
(gap  =  0.1  pm),  (c)  a  =  1.0  (gap  =  0.2  pm),  (d)  a  =  1.5  (gap  =  0.3  pm). 

HG.  2-  Exposure  trees  vs.  a  (dimensionless  gap)  for  Case  1.  (Gaps  range  from 
0-0.3  pm.)  The  line  is  biased  33%  (40  nm  resist  line),  (a)  MoM 
solution,  (b)  Rayleigh-Sommerfeld-Kirchhoff  approximation. 

FIG.  3  Intensity  vs.  spatial  frequency  for  Case  1.  Note  roll-off  at  0.05  nm*L 

FIG.  4  Intensity  vs.  lateral  position  for  Case  2.  (a)  a  =  0  (gap  =  0),  (b)  a  =  0.5 
(gap  =  3.7  pm),  (c)  a  =  1.0  (gap  =  7.5  pm),  (d)  a  =  1.5  (gap  =  11.2  pm). 

FIG.  5  Exposure  trees  vs.  a  (dimensionless  gap)  for  Case  2.  (Gaps  range  from 
0-11.2  pm.)  The  line  is  biased  50%  (150  nm  resist  line),  (a)  MoM 
solution,  (b)  Rayleigh-Sommerfeld-Kirchhoff  approximation. 

FIG.  6  Intensity  vs.  spatial  frequency  for  Case  2.  Note  roll-off  at  0.05  nm‘L 
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TABLE  I.  Computational  cases  coiisidered. 


Wavelength  X  (nm) 
Refractive  Index  [25] 
Width  W  (nm) 
Height  (nm) 
Transmission 
Phase  Shift  (rad) 
Patch  Size  (nm) 


Case  1 

4.48 

1  -  7.54x10-3  +  jl.04xlO-2 
30 
100 
0.0541 
1.058 
-0.64 


Case  2 

1.334 

1  -  2.31x10-3  +  jl.l9xl0-3 

100 

250 

0.0607 

2.72 

-0.19 
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3.1  Introduction 

For  microwave  integrated  circuit  applications,  the  characteristics  of 
intercoimects  have  been  investigated  for  propagation  modes  [1,2],  time 
response  [3],  crosstalk  [4],  coupling  [5],  delay  [6],  etc.  In  these  analyses, 
it  is  assumed  that  quasi-TEM  modes  are  guided  along  the  multicon¬ 
ductor  transmission  line.  In  [l],  the  analysis  was  performed  for  two 
asymmetric  transmission  lines.  In  [2]  and  [3],  an  arbitrary  number  of 
transmission  lines  were  analyzed.  In  [3],  the  load  and  the  source  condi¬ 
tions  were  presented  in  terms  of  the  modal  reflection  and  transmission 
coefficient  matrices. 

To  perform  the  quasi-TEM  analysis,  the  capacitance  matrix  for 
the  multiconductor  transmission  line  has  to  be  obtained  first.  Both 
the  spectral  and  the  spatial  domain  methods  have  been  proposed  to 
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S.  Modelling  of  Lossy  Microstrip  Lines 


Figure  8.1  Geometrical  configuration  of  M  microstrlp  lines  with  finite 
thickness  embedded  in  arbitrary  layers  of  an  isotropic  stratified  medium. 

calculate  the  capacitance  matrix.  In  the  spectral  domain  methods,  two 
side  walls  are  used  to  enclose  the  whole  transmission  line  structure, 
and  the  thickness  of  the  strip  lines  has  not  been  considered  [7,8],  In 
using  the  spatial  domain  method  [9],  the  structure  has  to  be  truncated 
to  a  finite  extent  to  make  the  numericed  implementation  feasible.  In 
[lO],  the  infinite  extent  of  the  structure  was  incorporated,  but  only  a 
two-layer  medium  was  considered. 

In  practical  microwave  integrated  circuits,  the  dielectric  loss  due 
to  the  substrate  and  the  conductor  loss  due  to  the  metallic  strips  have 
been  studied  in  the  analysis  of  circuit  performances[ll-13]. 

In  this  chapter,  we  present  a  quasi- TEM  analysis  of  coupled  lossy 
microstrip  lines  of  finite  strip  thickness  embedded  in  different  layers 
of  a  lossy  isotropic  stratified  medium  as  shown  in  Fig.  3.1.  First,  a 
spectral  domain  scalar  Green’s  function  in  a  lossy  isotropic  stratified 
medium  is  derived.  In  the  formulation,  no  side  walls  are  introduced, 
the  transmission  structure  is  not  truncated,  and  the  analysis  is  valid 
for  arbitrary  nvimber  of  dielectric  layers. 


3.2  Integral  Equation  Formulation 
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Based  on  the  scalar  Green’s  function,  a  set  of  coupled  integral 
equations  is  obtained  for  the  charge  distribution  on  the  strip  surfaces. 
The  method  of  moments  is  then  applied  where  pulse  basis  functions 
and  a  point-matching  scheme  is  used  to  solve  numerically  the  set  of 
integral  equations  for  the  charge  distribution,  and  hence  the  capaci¬ 
tance  matrix.  The  duedity  between  the  electrostatic  formulation  eind 
the  magnetostatic  one  is  applied  to  calculate  the  inductance  matrix. 
The  conductance  matrix  is  obtained  by  using  the  duality  between  the 
electrostatic  problem  and  the  current  field  problem.  A  perturbation 
method  is  used  to  Ccdculate  the  resistance  matrix. 

Finally,  a  transmission  line  analysis  is  derived  to  obtain  the  trans¬ 
fer  matrix  for  multiconductor  line,  which  significantly  reduces  the  effort 
in  treating  the  load  and  the  source  conditions.  Transient  responses  are 
obtained  by  using  the  Foxirier  transform.  The  results  for  two  coupled 
lines  are  presented. 


3.2  Integral  Equation  Formulation 


We  first  formulate  the  scaleir  Green’s  function  in  the  spectral  do¬ 
main  in  a  homogeneous  medium  of  permittivity  e  with  a  ground  plane 
located  at  z  =  0 .  Thus,  we  consider  a  uniform  line  cheirge  of  unit  am¬ 
plitude  to  be  located  at  along  the  y  direction,  and  evaluate 

the  electrostatic  potential  at  {x,z) . 

The  scalar  Green’s  function  ff(r,r')  is  the  solution  of  the  following 
Poisson  equation  : 

-  r')  (1) 

where  Vj  =  d'^ldz^  d^Jdz^  ,  f  =  iz  +  zz ,  and  f'  =  zx'  +  zz' .  The 
scalar  Green’s  function  can  be  expressed  in  the  spectral  domaun  as 


where  k,  =  rib,  -|-  zib,  ,  and  g{k,)  is  the  Fourier  transform  of  g{r,r') 
with  respect  to  x  and  z  .  Substituting  (2)  into  (1),  and  using  the  image 
theory,  we  get 


g{k,)  = 


47r^€(fc*  +  kl) 


— tfc,*'— tfc,  *' 


-*fc,  x'  +  i*,  i' 


(3) 


S.  Modelling  of  Lossy  Microstrip  Lines 


Figure  S.2  Definition  of  the  local  reflection  and  transmission  coefflcients 
at  a  dielectric  interface  between  media  (1)  and  (/  —  1). 


The  second  term  is  the  Fourier  transform  of  the  Green’s  function  due 
to  the  image  line  charge  at  f'j  =  xx'  -  zz'  which  is  the  image  of  r' 
with  respect  to  the  ground  plane  at  2  =  0.  The  scalar  Green’s  function 
can  then  be  written  as 


By  contour  integration  in  the  complex  ife,  plane,  we  have 


(4) 

(5) 


The  integrand  is  regular  at  =  0  ,  and  the  integration  is  well  defined 
as  r  goes  to  infinity. 

Next,  we  introduce  local  reflection  and  transmission  coefficients 
at  a  boundary  between  two  dielectrics.  Consider  a  dielectric  interface 
between  media  (f)  and  (I  -  1)  as  shown  in  Fig.  3.2.  We  assume  that 
the  spectral  domain  potential  in  each  region  can  be  expressed  as 


3.2  Integral  Equation  Formulation 


Figure  3.3  Definition  of  the  global  reflection  coefficients  at  the  upper 
boundaries  of  dielectric  layers  (/)  and  (I  —  1). 

where  zt  =  z/_i  =  z  +  di-i  ,  and  are  constants  to  be 

determined.  The  z  component  of  the  electric  field  in  each  region  can 
be  obtained  as 


Eu{f,  fc.)  [c-l''*'*'  -  (7a) 

(76) 

By  imposing  the  boundary  conditions  that  the  potential  and  the  nor¬ 
mal  component  of  the  electric  flux  density  are  continuous  at  z  = 
-dj_i ,  we  obtain 


g<  -  Cf-l 
ff  +  €<-l 
2et 

<{  +  ei-i 


(8a) 

(86) 


where  Ri(i-i)  and  can  be  interpreted  as  local  reflection  and 

transmission  coefficients,  respectively. 

Next,  consider  the  two  so\irce-free  layers  ( I )  and  ( /  -  1  )  of  the 
stratified  medium  as  shown  in  Fig.  3.3,  where  the  potential  and  the  z 
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3.  Modelling  of  Lossy  Microstrip  Lines 


component  of  the  electric  field  in  each  layer  can  be  expressed  as 

+  B,el''*l*']  (9a) 

Ej,{r,k,)  =£•''•*  (96) 

where  Zj  =  z  +  dj  and  j  =  / ,  /  -  1 .  We  introduce  the  global  re¬ 
flection  coefficient  defined  at  the  upper  boimdauy  of  layer  (j), 

z  =  -dj_i ,  as 

where  hj  is  the  thickness  of  layer  {j).  By  imposing  the  boundary 
conditions  at  z  =  -dj_i ,  we  obtain  the  following  recursive  relation 


1  =  2, .-.^N  (11) 


and  iZyi  =  iZio.  Similarly,  the  global  reflection  coefficient  Rnj  is 
defined  at  the  lower  boundary  of  layer  (j),  z  =  —dj  ,  as 

Rm  =  (12) 


By  imposing  the  boundary  condition  at  z  =  —  dj ,  we  obtain 


Rni  — 


Ri{l+\)  +  Rn(l+i)€  2|fc.l'‘i-n 

1  -I-  ’ 


/  =  iV-l,.--,0  (13) 


and  Rntf  =  —  1 . 

With  the  use  of  the  above  reflection  coefficients,  the  expression  of 
the  scalar  Green’s  function  in  a  layered  medium  can  be  obtained  in  a 
simple  way.  For  the  case  where  both  the  line  cheirge  and  the  observation 
point  are  located  in  the  upper  half-space,  that  is  layer  (0) ,  we  have 

1  /■'*>  e*k.  (*-»') 

|t.l 

(14) 


3.2  Integral  Equation  Formulation 


01 


where  zq  =  z  +  do,  rj  =  z'  +  do ,  the  uniform  line  charge  of  unit 
amplitude  is  located  at  r'  =  (z\zq),  and  the  observation  point  is 
located  at  f  =  (i,2o)  •  The  first  term  inside  the  brackets  is  the  direct 
term,  and  the  second  term  can  be  interpreted  as  due  to  the  reflection 
from  the  lower  bovmdary  at  z  =  —do  . 

For  the  case  when  both  the  line  charge  and  the  observation  point 
are  both  located  in  an  arbitrary  layer  (/)  with  1^0,  the  Green’s 
function  can  be  expressed  as 


(15) 


where  zi  =  z  +  di ,  z'l  =  z'  +  di ,  A{kg)  and  B{kg)  are  the  unknowns 
to  be  determined.  From  the  definition  of  Rui  and  Rr\i ,  we  have 


+  =  hi  (16a) 

iZn/ =  yl(Jfc,)e-l**l*',  r/  =  0  (166) 

By  solving  (16),  we  get 


1  - 


(17a) 

(176) 


Substituting  (17)  into  (15),  we  obtrdn  the  explicit  form  of  the 
Green’s  function  gu{r,T')  in  the  layer  (/)  as 


1 


+  iZn/e-l** ^ 


Figure  3.4  Determination  of  the  scalar  Green’s  function  where  the  source 
is  in  a  layer  (m)  and  the  observation  point  is  in  the  layer  (/).  (a)  Layer 
(/)  is  above  layer  (m).  (b)  Layer  (1)  is  below  layer  (m). 


The  first  term  inside  the  brace  is  the  direct  term,  the  other  terms 
are  the  summation  of  the  multiple  reflections  between  the  two  bound¬ 
aries  at  z  =  -dj_i  and  z  =  —dt .  It  is  observed  that  (14)  is  a  special 
case  of  (18)  with  1  =  0  and  R^i  =  0  . 

Next,  we  consider  the  case  when  the  source  and  the  observation 
points  are  in  different  layers.  As  shown  in  Fig.  3.4(a),  the  source  is 
in  a  layer  (m)  and  the  observation  point  is  in  a  layer  (/)  which  is 
above  layer  (m) .  We  assume  that  each  spectral  field  component  is 
transmitted  upward  from  the  source  layer  (m)  to  the  layer  (/)  with 
the  upward  transmission  coefficient  Xul,m  •  The  Green’s  fimction  can 
thus  be  expressed  as 

e|fc.lC  + 

1  - 


.’i* 


^  v;  -  *'*>  ■- 
.'*.,•  .-r  T-^^  '-'^- 
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where  zj  and  are  the  local  coordinates  defined  as  zj  =  z  +  dj  and 
=  z'  +  dm-i ,  respectively.  Similarly,  when  the  observation  point  is 
located  in  the  layer  (/  +  1) ,  we  have 

1  /■** 

«(w.)™(F>.f'<)  =  y__ 

cl**l*-  + 

where  Zj+i  =  z  +  dt+i  .  Imposing  the  boundary  conditions  at  the  in¬ 
terface  z  =  —  dj ,  we  get 

y  _  v  l^i+i  ^  •^*-<(1+1)  I  _  ~  o  n 


and  for  /  =  m  —  1 ,  we  have 


1  +  -Rum 


For  the  case  when  the  layer  (m)  is  above  layer  (/)  as  shown  in 
Fig.  3.4(b),  the  Green’s  function  may  be  expressed  as 

e-|fc»UU  + 

1  -  iZum.RnmC"^l*'*l'*"* 

where  Xnt,m  is  the  downward  trsmsmission  coefficient,  zj  and  are 
the  local  coordinates  defined  as  zi  =  z  +  dj_i  and  z^  =  z'  -f-  d^ , 
respectively.  Similarly,  when  the  observation  point  is  located  in  the 
layer  (/  -  1) ,  we  have 


I  fOO  (*-*') 

*  "  4^  7-00 

e“I**l*"  -f- 
1  - 


r-- ^  ^ 


>  »j*i‘>*.  *'*5  ^^.'.-•/  • 
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where  zj_i  =  z  +  (f/_2  .  Imposing  the  boundary  conditions  at  the  in¬ 
terface  z  =  —di-i ,  we  have 


Y  ,  —  X  /I  ■,\  (;~l**«l**‘-»  ^  ^n(t-i) 

Ani.m  -  An(j-i).me  ^  ’ 


and  for  /  =  m  d-  1 ,  we  have 


I  =  m  +  2,--.,iV 
(25) 


1  +  -fir 


Next,  consider  M  microstrip  lines  embedded  in  arbitrary  layers  of 
an  isotropic  stratified  medium  as  shown  in  Fig.  3.1.  All  the  conductor 
strips  are  of  finite  thickness,  the  permittivity  of  layer  (j)  is  ej  ,  and  a 
groimd  plane  is  put  at  z  =  ~djf  as  the  potentizJ  reference. 

Using  the  scalar  Green’s  function,  the  potential  in  layer  (/)  can 
be  represented  as 

M 

^K^)  =  X)/  ‘^9lL(p){ryr)pp{f')  (27) 

P=l 

where  Fp  is  the  cross  section  contour  of  the  pth  microstrip  line,  Pp{r*) 
is  the  charge  density  on  the  pth  microstrip  surface,  L{p)  is  the  layer 
where  the  p  th  microstrip  line  is  embedded.  The  cross  section  of  the 
microstrip  lines  can  be  arbitrary  in  genersd.  For  practiced  applications, 
ordy  the  rectangular  cross  sections  ewe  considered. 

The  potential  distribution  in  (27)  satisfies  all  the  boundary  condi¬ 
tions  that  the  potential  and  the  normal  electric  flux  density  are  con¬ 
tinuous  across  the  interfaces  between  adjacent  dielectric  layers.  To  ob- 
tedn  the  charge  distribution  on  each  microstrip  surface,  we  impose  the 
boimdary  condition  that  the  potential  on  each  microstrip  is  equal  to 
the  impressed  voltage.  Thus  we  have 


M  . 

^  /  dr' ffL(,)i,(p)(r,r')pp(f')  =  V„ 

p=i  ■'r. 


r  on  Fa 


where  g  =  •,M. 

In  the  next  section,  the  method  of  moments  is  applied  to  solve  (28) 
for  the  charge  distributions. 


S.S  Numerical  Solution  of  Charge  Distribution  Bo 

3.3  Numerical  Solution  of  Charge  Distribution 

By  applying  the  method  of  moments  to  solve  (28),  we  first  choose 
a  set  of  pulse  basis  functions  for  the  surface  charge  density,  hence 

=  p  =  l,-.-,M  (29) 

t=i 

where  Up,-  and  Ppiif')  are,  respectively,  the  expansion  coefficients  and 
the  basis  functions  on  the  ith  section  of  the  pth  microstrip  surface; 
and  Np  is  the  total  number  of  basis  functions  on  the  pth  microstrip 
surface.  The  pulse  function  is  defined  as 


1,  r  on  Fp,- 
0,  elsewhere 


where  Fp^  is  the  i  th  section  on  the  p  th  microstrip  surface. 
Substituting  (29)  into  (28),  we  have 


M  N, 


T  on  Fo 


iTS  •'W  m 

=  V,,  f 

p=i  «=i 


where  q  =  1 ,  •  •  • ,  M . 

Next,  we  choose  the  center  point  at  each  pulse  basis  as  the  testing 
point  to  test  (31),  thus  we  have 

M  Nf 

=  <iV,  (32) 

i  =  l 


where 


Jr  mi 


=  V,  (34) 

and  =  (z,j,2,j)  is  the  center  coordinate  of  F,j  .  By  solving  (32) 
for  Opi ’s,  the  charge  density  can  be  obtained. 

The  capacitance  matrix  U  of  M  microstrip  lines  can  be  obtained 
by  solving  (28)  M  times.  For  the  n  th  time,  all  the  microstrip  lines  are 
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grounded  except  that  one  volt  is  imposed  on  the  n  th  line.  By  solving 
(32),  the  charge  distribution  is  obtained.  The  total  surface  charge  per 
unit  length  on  the  m  th  microstrip  line  is  equal  to  the  mn  th  element  of 
the  capacitance  matrix,  Cmn  ■  This  can  be  observed  from  the  definition 
of  the  capacitance  matrix 

M 

=  l<g<M  (35) 

p=\ 

By  setting  Vp  =  5p_n ,  vire  have  with  I  <  q  <  M .  Here, 

Sp  n  is  the  Kronecker’s  delta  fimction  and  the  superscript  (n)  in 
is  the  index  of  the  microstrip  where  one  volt  is  imposed. 


3.4  Magnetostatic  Dual  Problem 


In  this  section,  we  briefly  review  the  magnetostatic  dual  problem. 
For  nonmagnetic  materials,  the  magnetostatic  problem  for  the  geomet¬ 
rical  configuration  shown  in  Fig.  3.1  is  equivalent  to  the  one  where  all 
the  dielectric  media  are  replaced  by  the  free  space.  The  magnetostatic 
potential  tp(f}  due  to  a  uniform  line  current  at  f  can  be  derived  in 
a  way  similar  to  that  in  Section  2  as 


p=i 


r. 


=  v-,,  r  on  r. 


(36) 


where  q  —  r')  is  the  scalar  Green’s  function; 

Jp(r')  is  the  surface  current  on  the  pth  microstrip  surface;  and  rpq  is 
the  magnetostatic  potential  measured  on  the  q  th  microstrip  surface 
with  respect  to  the  groxmd  plane,  which  is  equal  to  the  magnetic  flux 
linkage  between  the  q  th  microstrip  line  and  the  ground  plane  per  unit 
length. 

Instead  of  solving  (36)  directly,  a  dual  electrostatic  problem  for 
the  same  geometrical  configuration  in  free  space  is  solved.  It  can  be 
shown  that 

9'^ir,r)  =  f^o(o9^ir,r')  (37) 

where  is  the  scalar  Green’s  function  derived  in  Section  2  with 

all  £; ’s  replaced  by  the  dielectric  constant  in  the  free  space. 


f  ^ iA  ;  .’''<• 


«■  !*-.  •  Wi.  ,'.  -'■  * 
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Imposing  the  same  boundary  conditions  that  V,  =  ^,  for  I  <  q  < 
M ,  we  have 

7W(r)=^^,  l<p.jb<M  (38) 

Hofo 

where  Pp*^(r)  is  the  stirface  charge  distribution  on  the  p  th  microstrip 
surface  when  we  impose  one  volt  on  the  k  th  strip,  and  zero  volt  on 
the  other  strips;  jj!'\r)  is  the  surface  current  distribution  on  the  pth 
microstrip  surface  when  we  impose  one  tesla-m  on  the  k  th  strip,  and 
zero  tesla-m  on  the  other  strips. 

In  genera],  the  surface  current  on  the  pth  microstrip  surface  can 
be  expressed  as 

M')  =  E  =  rV  L  '►‘Op ’(f )  (39) 

Integrating  (39)  over  one  unit  length  on  the  pth  microstrip  surface, 
we  have 

1  ^ 

/p = — y2 

where  Ip  is  the  surface  cxirrent  on  the  p  th  microstrip  line,  and  Co,pfc  is 

the  pk  th  element  of  the  capacitance  matrix  Co  with  all  the  dielectrics 
replaced  by  free  space.  Hence,  we  obtain  the  conventional  result  that 
the  inductance  matrix  L  is  proportional  to  the  inverse  of  the  capaci¬ 
tance  matrix  in  free  space  as 

L  =  po^o^o  (41) 


Calculation  of  ^  and  R 
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where  the  conductance  matrix  G  accounts  for  the  dielectric  loss,  and 
the  resistance  matrix  R  accounts  for  the  conductor  loss. 

The  nth  eigen  solution  to  (42)  can  be  represented  as 

7  =  7„e-''-«'  (43a) 

V  =  Vne""""  (436) 

where  7„  is  the  eigenvalue  of  the  nth  mode;  and  /„  and  Vn  are  the 
corresponding  eigenvectors.  Substituting  (43)  into  (42),  we  have 


7n/n  =  [c?  -  (44o) 

7„v„=  (446) 

and  hence 

7^7„  =  p  -  .  [1  -  twl]  .  (45a) 

7*Fn  =  [I  '  •  p  -  (456) 

There  are  in  general  M  eigen  solutions  to  the  above  equations.  The 
eigenvaJues  solved  from  (45a)  and  (45b)  are  the  same;  and  the  eigen¬ 
vectors  In ’s  and  Vn ’s  are  related  by  (44). 

For  the  n  th  eigenmode,  the  time-average  Poynting’s  power  P„r 
can  be  represented  as 

Pnr  =  ^Re[7l'V„]  (46) 

where  /],  is  the  transposed,  complex  conjugate  of  •  The  time- 
average  power  loss  per  unit  length  PnL  along  the  transmission  lines 
can  be  calculated  as 


PfiL  =  PnC  +  PnD  =  2a„rPnr  (47) 

where 

PnC  =  (l/2)7l  •  1  •  7„  =  2ancPnT  (48a) 

PnD  =  (l/2)Vl  ■^’Vn  =  2anDPnT  (486) 


■.r.  ■  .  i: 


a  •  .xS:-rZ  It 
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where  Pnc  P„£>  are  the  time-average  conductor  loss  and  dielec¬ 
tric  loss  per  unit  length,  respectively;  a„r ,  a„c .  and  a„£>  are  the 
corresponding  attenuation  constants  for  the  n  th  mode. 

To  calculate  the  conductance  matrix  which  accounts  for  the  dielec¬ 
tric  loss,  we  make  use  of  the  duality  between  the  electrostatic  problem 
and  the  current  field  problem.  Assume  that  the  geometrical  configu¬ 
ration  of  the  current  field  problem  is  the  same  as  in  Fig.  3.1,  and  the 
conductivity  in  layer  (1)  is  designated  by  <r/  with  0  <  /  <  iV  .  If  the 
conductivities  in  the  current  field  problem  and  the  dielectric  constants 
in  the  electrostatic  problem  satisfy  the  following  relation 


<ro  <Ti 


then,  we  have 


The  perturbation  method  is  used  to_solve  for  the  attenuation  con¬ 
stants  a„c  and  the  resistance  matrix  R .  We  start  from  the  lossless 
eigenvalue  equations  : 


C  '  L  •  lon  =  —^lon 
vl 


T  •  C  •  Von  =  “5"  I' On 


where  7„  =  -i0n  ;  Vn  =  ^IPn  is  th®  phase  velocity  for  the  n  th  mode; 
and  C  is  the  capacitance  matrix  with  all  materizds  lossfree.  The  time- 
average  guided  power  can  be  approximated  by  the  power  guided  sJong 
lossless  lines  as 

Pnr«  (52) 

The  power  loss  per  unit  length  of  the  nth  mode  due  to  imperfect 
conductor  can  be  calculated  by  using  the  surface  current  obtained  for 
the  perfect  conductor  lines.  Thus,  we  have 


where  Rks  is  the  surface  resistance  per  unit  area  on  the  k  th  mi¬ 
crostrip  surface,  Jnh{r)  is  the  surface  current  distribution  of  the  nth 


- 

'-'.i-C' “l'  n‘ 


•  ‘  ®  I 

■M  ■■ 
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mode  on  the  k  th  microstrip  siirface.  Here  M'  is  the  total  number  of 
conductors  including  the  grovmd  planes.  As  defined  in  (29),  cam 

be  expanded  in  terms  of  the  pulse  basis  functions  as 

p'llhf)  =  E  (54) 

m=l 

where  Pkni^)  is  the  basis  function  defined  in  (30).  By  using  (39),  the 
surface  current  on  the  k  th  microstrip  surface  can  be  represented  as 

•'‘(f)  =  rV  S  *  E  42.r*™(f )  (55) 

w‘«  fe; 

where  tpi  is  obtadned  by  solving 

M 

V**  —  ^  ^  (^®) 

i=i 

where  Lij  is  the  ij  th  element  of  the  inductance  matrix  L ,  Ion,j  is 
the  j  th  element  of  the  eigenvector  lon  •  The  conductor  loss  can  then 
be  calculated  as 

.  M  Nt  M  *  . 

^nc  =  2  E  ^<«5(/^OCo)-*  E  E  /  inm(r)|*  dr  +  Po 

fc=l  m=l  i=l  •'I'*"* 

(57) 

where  Pc  is  the  conductor  loss  due  to  the  grovmd  plames.  To  calculate 
Pa  ,  we  first  solve  for  the  magnetic  field  Hg  as 

fc=im=iLi=i  J  <'r»« 

7-00  ^  Po  dz  J 


where  p^(Jfc*,  2,  z')  is  the  Fourier  transform  with  respect  to  *  of 
p^(r,r')  as  defined  in  (37).  Thus,  we  have 

M  Ni  M  T  M  1  •  r  M 

Pa  =  2tRs  E  £  E  E  E 

fc=l  m=l  V=1  m'=l  L»=l  L«=l 


,  f  c"  --iVV,.  -A.  V  -•  >  >-.*.-:t-5^4fcv?)>| 


■<i.YS 


•'  V  -i,  •. 
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f  dfPkmir')  f  dr"Pk-rn‘{^')  f  dfc,  cos[fc.(z' -  x")] 

r  1  ^G» ^*)l  *  r _ 1  dg^{kxt2Q,z  )i  ,gg. 

I  Un  dz  J  I  lift  dz  J 


where  Rs  is  the  surface  resistance  per  unit  area  of  the  ground  planes, 
zq  is  the  z  coordinate  of  the  ground  planes. 

Now,  set  G  =  0  in  (45b),  we  have 


(a„c  -  iPn?Vn  =  [ii  -  •  [-iwC]  •  V 


Making  the  approximation  that  /3„  a  /?0n »  V^n  **  ^On »  utilizing 
the  relation  (44)  with  G  =  0 ,  the  imaginary  part  of  (60)  becomes 


2£rnC^n,m  —  ^  ^  ■Rmfc-fpft.fct  1  —  Tl,  TTl  <  M 


where  /on.fc  i*  the  k  th  element  of  lon »  and  Vbn.m  is  the  m  th  element 
of  Von  •  The  elements  of  R  can  1«  oMai^d  by  solving  (61). 

After  obteiining  the  matrices  ^ ,  L  ,  G ,  and  R  ;  the  matrix  equa¬ 
tions  (44)  and  (45)  can  be  solved  for  the  eigenveJues  and  eigenvectors. 


3.0  Ttansmission  Line  Analysis 


The  voltage  and  current  along  the  transmission  lines  can  be  rep¬ 
resented  in  terms  of  the  eigenvectors  obt^ned  in  the  last  section  as 


F(y)  =  a„Vn8-^-‘'  +  £  f>nV„c-' 


=  Sv  ’  A(y)  •  A  +  Syr  •  A(— y)  •  B 


(62a) 


/(»)  =  E  -  E  W"'" 


=  Sj  •  A(y)  •  A  —  Sj  •  A(— y)  •  B 
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where  5v  is  the  matrix  with  the  eigenvector  Vn  as  its  nth  column; 
Sj  is  the  matrix  with  the  eigenvector  as  its  n  th  column;  A  and 
B  axe  the  mode  coefficient  vectors  in  the  forward  and  backward  direc¬ 
tions,  respectively;  A(y)  is  a  diagonal  matrix  with  exp(— 7„y)  as  its 


n  th  element;  namely, 

"X  =  (oi,  oj,  •  •  • ,  (63a) 

_B  =  {bub2,^-‘,bMf  (636) 

A(y)  =  diag.  •  •  • ,  (63c) 

where  the  superscript  T  represents  the  transpose  of  a  row  vector. 

The  characteristic  impedance  matrix  can  be  defined  as 

Zc  —  Sv  "  Sj  (64) 

From  (62),  the  mode  coefficients  A  and  B  can  be  represented  in 
terms  of  the  line  voltages  V(0)  and  line  currents  7(0)  as 

^=5pv'*V"(0)  +  f7'*7(0)]  (65a) 

^pv'-V^(0)-fj''7(0)]  (656) 


Next,  the  line  voltages  y(/)  and  the  line  currents  /(/)  can  be 
expressed  in  terms  of  y(0)  and  7(0)  through  (65)  as 

no]  r3vv(/)  iv/(/)]  rmi  _  moy 

=  _  _  •  _  =^(0-  _  (66) 

[ami)  Au{i)j  lj'co)]  Ino), 

where  A{I)  is  called  the  treinsfer  matrix  for  the  uniform  treinsmission 
lines  of  length  /;  and  the  explicit  form  of  the  submatrices  of  A(I)  are 


^vv(0  =  •  [X(0  +  ^(-o]  • 

(67a) 

^v/(o  =  ^^K*p(o-^(-o]-^7' 

(676) 

%v{l)=pi’^{l)-A{-l)]-Ty 

(67c) 

I//(0  =  5^/-p(0  +  ^(-o]-^/"' 

(67d) 

3.0  Tranamiuion  Line  Analysis 
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If  more  than  one  transmission  line  section  with  different  characteristics 
are  cascaded,  the  tramsfer  matrix  of  each  section  is  multiplied  to  obtain 
the  overall  transfer  matrix.  _ 

Imposing  a  voltage  source  Vy  in  series  with  an  impedance  matrix 
Zy  at  y  =  0,  and  imposing  a  voltage  source  Vl  iu  series  with  an 
impedance  matrix  Zl  at  y  =  / ,  we  have 


r(0)  =  Fy-Zy7{0) 

v{i)  =  m 

The  line  voltages  at  all  ports  can  thus  be  calculated  as 

7  Aviil)'Zs  -  Avv(l) 


(68a) 

(686) 


Zr 


Avi{i)-Zs  '^s 


7(0). 

An{l)-Zs  'Vs  +  Zi  -Vl,} 


(69) 


The  transient  response  can  be  obtained  by  using  Fourier  transform. 

For  the  case  of  a  single  microstrip  line,  we  first  solve  for  the  ca¬ 
pacitance  and  the  inductance  per  unit  length  C  and  L ,  respectively. 
Next,  solving  the  eigenvalue  equations  in  the  lossless  medium,  we  have 


7*  =  =  -uHC 

/o  =  l 


Vo=  y/L/C 


(70a) 

(706) 

(70c) 


The  magnetic  flux  linkage  is  ^  =  LIq  =  L ,  the  time-average  power 
is  Px  =  (\l2)y/LIC ;  and  the  conductor  loss  per  unit  length  can  be 
calculated  by  (57)  and  (59),  hence  we  have  ac  =  PcI“^Pt-  Using 
(61),  we  have  R  =  2ac7  y/LjC .  The  conductance  per  unit  length  can 
be  calculated  by  using  (50).  With  L,  C ,  R,  and  G ,  the  eigenvalue 
and  the  eigenvector  can  be  solved  from  (44)  and  (45)  as 


7  =  V(G  -  iu)C){R  -  ia>L) 
1  =  1 


(71a) 

(716) 


mE:;^ 


104  S.  Modelling  of  Lossy  Microstrip  Lines 

The  transfer  matrix  A{1)  for  a  transmission  line  of  length  /  can  be 
obtained  from  (67)  as 


^(/)  = 


cosh?/  -Z£7sinh7/ 
—2c^siiihyl  cosh?/ 


For  the  case  of  two  symmetric  microstrip  lines,  we  first  solve  the 
eigenvalue  equations  in  the  lossless  medium  for  the  even  and  the  odd 
modes,  respectively.  Thus  we  have 

?e  =  —fie  =  4-  iij)(Cii  +  Cij)  (73a) 


Ie0  = 


17  _  /  T'li  +  -^17 
1  ’  + 


(73a) 


7o  —  '  ■^12)(C'll  ~  Cli) 

/.«=  f  1 


I  Lii  —  Li2 

('ll  -  Cl2 


^  L\\  —  Xia  j 

Cn-C„  “ 


(73c/) 


The  magnetic  flux  linkage  and  the  time- average  power  tan  be  calcu¬ 
lated  as 


xp^  =  L  ‘  leo  =  (^11  +  Li2)Ie0 

(74a) 

„  l^ll  +  Li2 

"  V  Cu  -1-  C22 

'746) 

xpg  —  L-  loo  —  (f-ii  —  ^12)  -ToO 

(74c) 

„  /in -1-12 

V  c„  -  C„ 

(74d) 

The  conductor  loss  per  unit  length  can  be  obtained  by  using  (57) 
and  (59),  thus  we  have 

QieC  =  i^«c/2P«T.  OLoC  =  PoCl’^PoT  (75) 

The  resistance  matrix  can  be  solved  from  (61)  as 


^11  =  a«c;i 
7^12  =  <^tC  ‘ 


L\i  -1-  L\2  . 

/ -^-ll  -  I'\2 

c-„  +  c.,  +  “'^\ 

'  Cii  -  Ci2 

Lii  L\2 

iLii  -  L\2 

Cii  -1-  C12 

1  Cu  -  C\2 

(76a) 

(76b) 
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The  conductance  matrix  per  unit  length  can  be  cjilculated  by  using 
(50).  Substituting  these  C,  R,  and  G  into  (44)  and  (45),  we 
obtain  the  eigenvalues  and  eigenvectors  for  the  even  and  the  odd  modes 
as 


7e  =  (yii  +  yi2)(2ii  +  ^ij) 


/e  = 


V,  =  Ze/. 


il  =  (yii  -  »i2)(^u  -  212) 


Io  = 


1 

-1 


Vo  =  Zolo 


(77a) 

(776) 

(77c) 

(77d) 


where 


z.  = 


y\\  +  yi2 


-  fn 
yi2 


(78) 


i'll  =  G\i  —  tU/Cii,  1/12  =  G\2  ~  i^Cl2 

2=11  =  Rn  —  iuLiif  Z12  —  R\2  ~ 


(79) 


Zc  =  Sv  •  Sj  = 


(80) 


The  characteristic  impedance  matrix  can  be  obtained  from  (64)  as 

(Z. +  Z„)/2  (Z,-Z„)/2' 

(Ze-Zo)/2  (Z,  +  Z„)/2 

The  transfer  matrix  of  length  I  can  be  calculated  from  (67)  as 
cosh7el  +  cosh7ol  cosh7el  —  cosh7ol 
cosh  7e^  —  cosh  7o7  cosh  7*/  +  cosh  7^/  J 


Avv{l}  =  2 


(81a) 


Aviil)  =  \ 

—  Zf  sinh  7*1  -  Z^  sinh  7^/ 
—  Z«  sinh  7,1+  Z^  sinh  70/ 


—  Ze  sinh7e/  +  Z^  sinh  70I 

—  Z,  sinh  7,1  —  Zo  sinh  7©! 


(816) 
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Aviil)  =  - 


-Z~^  sinh7e/  -  Z~^  sinh^o^  -^e  *  sinh7e/  +  Z~^  sinh7^/] 

sinh7,/  +  Z~^  sinh7of  -Z~^  sinh7el  -  Z~^  sinh7o/ 

(81c) 


^Mi)  =  j 


cosh7*/ +  cosh7ol  cosh  7*  1  -  cosh  70/ 
cosh  7e/ —  cosh  7o/  cosh  7*  f  +  cosh  7,,/ 


(81d) 


In  the  next  section,  we  demonstrate  how  the  transfer  matrix  is  applied 
to  solve  the  coupled  transmission  line  problems. 


3.7  Results  and  Discussions 

First,  we  check  the  results  of  our  method  in  calculating  the  capac' 
itance  matrix  using  the  spectral  domain  Green’s  function  and  compare 
it  with  other  methods  using  the  spatial  domain  Green’s  function.  In 
Tables  3.1  to  3.4,  we  present  the  results  of  capacitance  matrix  for  differ¬ 
ent  microstrip  line  configurations.  In  Table  3.1,  the  difference  between 
our  results  and  those  in  [9]  is  about  1%  for  the  self-capacitance,  and 
is  about  0.03%  for  the  mutual  capacitance.  The  difference  from  those 
in  [10]  is  about  0.1%  for  the  self-capacitance,  and  about  1%  for  the 
mutual  capacitance. 

In  Table  3.2,  the  difference  between  our  results  and  those  in  [9]  is 
less  than  0.7%  for  the  self-capacitance,  and  is  less  than  4%  for  the 
mutual  capacitance.  The  difference  between  our  results  and  those  in 
[10]  is  about  0.01%  for  the  self-capacitance,  imd  is  less  than  0.4%  for 
the  mutual  capacitance. 

The  capacitance  matrix  for  two  microstrip  lines  embedded  in  the 
saune  and  different  layers  of  a  two-layered  medium  are  presented  in 
Tables  3.3  and  4,  respectively.  The  difference  between  our  results  and 
those  in  [9]  is  about  3%  for  the  self-capacitance  and  about  1%  for  the 
mutual  capacitance. 

In  Fig.  3.5,  we  present  the  charge  distributions  on  the  broad  sides 
of  a  microstrip  line  with  t/h  =  0.02  and  w/h  =  0.1 ,  1.0,  2.0.  We  use 
24  pulses  on  each  broad  side,  and  2  pulses  on  each  narrow  side.  It  is 
observed  that  the  charge  density  on  the  bottom  side  of  the  microstrip 


3.7  Results  and  Discussions 


CoznpsriioD  of  Capsdtsnce  Matrix  Elemeots  (F/m). 


clement 

present  work 

(«] 

|10) 

0.6301  X  10-** 

0.6233  X  10“'* 

1  0.6307  X  10"'* 

^15 

-0.5929  X  10-” 

-0.5931  X  10“" 

-0.5866  X  10"" 

Cj2 

0.6301  X  10-'* 

0.6233  X  10-'* 

0.6307  X  10"'* 

Table  3.1  Comparison  of  capacitance  matrix  for  two  symmetrical  strip* 
lines  of  Unite  thickness. 


is  about  an  order  larger  thain  that  on  the  top  side  because  the  electric 
field  between  the  microstrip  and  the  ground  plane  is  stronger  than  the 
electric  field  above  the  microstrip.  Also,  as  the  width  of  the  microstrip 
increases,  the  edge  effect  becomes  less  significant. 

In  Table  3.5,  we  present  the  results  of  resistance  calculation  for  a 
microstrip  line  compared  with  those  in  [12].  The  discrepancy  is  7.6% 
for  w/h  =  0.1,  and  3.4%  for  w/h  —  2.0.  The  calculation  of  the 
resistance  depends  on  the  square  of  the  charge  density  distribution 
which  possesses  edge  effect  as  shown  in  Fig.  3.5.  Hence,  even  when  two 
different  methods  can  predict  close  capacitance  results,  it  is  possible 
that  the  resistance  results  can  deviate  by  a  higher  percentage. 

Next,  we  present  the  frequency  response  and  the  transient  response 


'  tz  I 


-f  '1 

■■ 
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elemtnt 


Compsrisoo  of  Capacitance  Matrix  Elements  (F/m). 


[10] 


Table  3.2  Comparison  of  capacitance  matrix  for  two  symmetrical  mi¬ 
crostrip  lines  of  finite  thickness. 

of  two  symmetrical  microstrip  lines.  The  driving  voltage  is  assumed  to 
have  a  sinusoidal  pulse  waveform  Vi(t)  with  duration  r  =  200  pico 
seconds  as 


Vi{t)  = 


(1/2)[1  -  cos(27rt/T)],  0  <  t  <  r 

0,  elsewhere 


The  frequency  response  and  the  transient  response  are  presented 
in  Fig.  3.6(a)  and  Fig.  3.6(b),  respectively,  with  all  the  four  load 
impedances  equal  to  the  characteristic  impedance  of  the  even  mode 
ZeQ.  In  Fig.  3.6(a),  Vi(0)  is  the  voltage  at  z  =  0  along  line  1,  Vi{l) 
is  the  voltage  at  the  receiving  port,  V2(0)  is  the  near-end  coupling, 
and  V2{1)  is  the  far-end  coupling.  The  increase  of  V2(/)  with  fre¬ 
quency  shows  that  the  high  frequency  components  are  responsible  for 
the  far-end  coupling[l6]. 
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Comparison  of  Cspscilance  Matrix  Elements  (F/m). 


element 


present  work 

l»l 

0.3E27  X  10- “ 

0.3720  X  10-‘» 

-0.6884  X  10-*> 

-0.6889  X  10-“ 

0.2245  X  10-** 

0.2169  X  10-‘* 

-0.3  .0.1  0.1  0.3 


Table  3.3  Comparison  of  capacitance  matrix  for  two  microstrip  lines  of 
finite  thickness  embedded  in  the  same  layer  of  a  two-layered  medium. 


In  Fig.  3.6(b),  the  even  and  the  odd  modes  propagate  in  different 
velocities,  and  the  odd  mode  propagates  faster  than  the  even  mode. 
Hence,  the  waveform  Vi{l)  becomes  broader  than  Vi(0) ,  and  V2({) 
shows  the  split  of  the  odd  and  the  even  modes.  The  waveforms  V^i(O) 
and  F2(0)  after  t  =  800  ps  are  due  to  the  r  eflections  by  the  load 
imped2tnces  at  y  =  / .  Since  the  load  impedances  are  chosen  to  be  the 
same  as  the  characteristic  impedance  of  the  even  mode,  only  the  odd 
mode  is  reflected. 

In  Fig.  3.7,  the  transient  response  is  presented  with  all  the  four 


m 
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-0.1  0.1 


<T  — *  OO 

Table  3.4  Comparison  of  capacitance  matrix  for  two  microstrip  lines  of 
finite  thickness  embedded  In  different  layers  of  a  two-layered  medium. 

load  impedances  equzd  to  the  characteristic  impedance  of  the  odd  mode 
^oo  •  The  split  of  the  odd  and  the  ev^n  modes  is  observed  again.  The 
waveforms  yi(0)  and  ^2(0)  after  t  =  800  ps  show  that  only  the  even 
mode  is  reflected,  and  the  reflected  signal  arrives  later  than  that  in 
Fig.  3.6(b)  because  the  velocity  of  the  even  mode  is  slower  than  that 
of  the  odd  mode. 

In  Fig.  3.8  and  Fig.  3.9,  we  present  the  transient  responses  of 
two  symmetrical  microstrip  lines  with  3  =  0.25  mm  and  a  =  0.375 
mm,  respectively.  The  load  impedances  are  chosen  the  same  as  the 
characteristic  impedance  of  the  even  mode.  Compared  with  Fig.  3.6(b), 
it  is  observed  that  the  coupling  signals  ^2(0)  and  V2{1)  become  weaker 
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•0. 5-0. 4  -0. 3*0. 2*0.1  0.0  0.1  0.2  0.3  0.4  0.5 

xfw 

Figure  S.6  The  charge  distribution  on  the  broad  sides  of  a  microstrip 
line  with  finite  thickness,  i,  =  11.7,  h  —  2cm,  t/h  —  0.02,  potential  on  the 
microstrip  surface  la  1  volt,  24  pulses  per  broad  side,  and  2  pulses  per 
narrow  side. 

as  the  sepairation  s  is  increased. 

In  Fig.  3.10,  the  transient  responses  with  a  complex  dielectric  con¬ 
stant  c,  =  10  -f  tO.l  are  presented.  The  sign  d  amplitudes  are  smaller 
than  those  in  Fig.  3.6(b)  due  to  the  dielectric  loss.  In  Fig.  3.11,  the 
transient  responses  with  copper  as  the  conductor  material  are  pre¬ 
sented.  The  surface  resistance  is  assumed  to  be  2.61  x  10“^\/7  ohms. 
The  signal  amplitudes  are  slightly  different  from  those  with  a  perfect 
conductor  because  the  copper  itself  is  good  conductor. 
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Comparison  of  Conductor  Loss  Parameters. 

w/h 

R(i\/m) 

^(n/m)(l2) 

aZoK/R,  (dB) 

aZoh/R.  (dB)[12] 

0.1 

0.07125 

0.06605 

23.713 

21.981 

0.2 

0.04477 

0.03993 

14.899 

13.289 

0.3 

0.03398 

0.02986 

11.309 

9.937 

0.4 

0.02795 

0.02400 

9.303 

8.120 

0.6 

0.02124 

0.01848 

7.069 

6.150 

1.0 

0.01492 

0.01315 

4.967 

4.376 

1.2 

0.01308 

0.01173 

4.353 

3.904 

1.4 

0.01165 

0.01063 

3.876 

3.538 

2.0 

0.00872 

0.00843 

2.901 

2.807 

<0 


.w. 


<r€o 


a  —*  ao 


Table  3.5.  Comparison  of  the  conductor  loss  parameters  for  a  microstrip 
line  of  finite  thickness. 


Conclusions 

The  spectral  domain  scalar  Green’s  function  in  a  lossy  isotropic 
stratified  medium  is  derived.  A  rigorous  integral  equation  formulation 
for  the  charge  distribution  on  the  surfaces  of  the  microstrip  lines  with 
finite  thickness  embedded  in  arbitrary  layers  of  an  isotropic  stratified 
medium  is  derived.  Using  the  spectral  domain  Green’s  function,  a  mul¬ 
ticonductor  transmission  line  analysis  is  formulated  to  investigate  the 
propagation  properties  of  coupled  lossy  microstrip  lines.  Both  the  fre¬ 
quency  and  the  transient  responses  of  coupled  lines  with  different  load 
conditions  can  be  obtained.  An  efficient  algorithm  is  devised  based  on 
this  approach. 
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Figure  3.8  Transient  response  of  two  symmetric  microstrip  lines  h  = 
0.2mm,  w  =  0.125mm,  t  —  5pm,  s  —  0.25mm,  I  =  5cm,  =  10,  all 
conductors  are  perfect,  Zi  =  Zt  ^  Zt  =  Z,o  ^  66.3450. 


a  268  460  608  660  1660 
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Figure  3.9  Transient  response  of  two  symmetric  microstrip  lines,  h  — 
0.2mm,  w  =  0.125mm,  t  =  5pm,  s  =  0.375mm,  I  =  5cm,  e,  =  10,  al 
conductors  are  perfect,  Z\  =  Zi  Z%  =  Z^  —  Z,o  %  63.5670. 
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Figure  S.IO  Transient  response  of  two  symmetric  microstrip  lines,  h  = 
0.2mm,  w  =:  0.125mm,  t  =  5pm,  $  =  0.125mm,  I  =  5cm,  10  +  >0.1,  all 
conductors  are  perfect,  Zi  =  Zi  =  Zt  =  Z4,  =  Z,o  ^  72.2700. 
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Figure  S.ll  Transient  response  of  two  symmetric  microstrip  lines,  h  = 
0.2Run,  to  =  0.125mm,  t  =  5pm,  «  =  0.125mm,  /  =  5cm,  (r  =  10,  al 
conductor*  are  copper,  Zi  =  Zi  =  Zt  =  Z^  =  Z,o  %  72.2700. 
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A  Hybrid  Method  for  the  Calculation  of  the 
Resistance  and  Inductance  of  Transmission 
Lines  with  Arbitrary  Cross  Sections 

Michael  J.'Tsuk,  Member,  IEEE,  and  Jin  Au  Kong,  Fellow,  IEEE 


Abstract — The  frequency-dependent  resisUnce  and  induc¬ 
tance  of  uniform  transmission  lines  are  calculated  with  a  hybrid 
technique  that  combines  a  cross-section  coupled  circuit  method 
with  a  surface  integral  equation  approach.  The  coupled  circuit 
approach  is  most  applicabie  for  low-frequency  calculations,  while 
the  integral  equation  s  pproach  is  best  for  high  frequencies.  The 
low-frequency  method  consists  in  subdividing  the  cross  section 
of  each  conductor  into  triangular  filaments,  each  with  an  as¬ 
sumed  uniform  current  distribution.  The  resistance  and  mutual 
inductance  between  the  filaments  are  calculated,  and  a  matrix  is 
inverted  to  give  the  overall  resistance  and  {t<r<actance  of  the 
conductors.  The  high-frequency  method  expresses  the  resistance 
and  inductance  of  each  conductor  in  terms  of  the  current  at  the 
surface  of  that  conductor  and  the  derivative  of  that  current 
normal  to  the  surface.  A  coupled  integral  equation  is  then 
derived  to  relate  these  quantities  through  the  diffusion  equation 
inside  the  conductors  and  Laplace’s  equation  outside.  The 
method  of  moments  with  pulse  basis  functions  is  used  to  solve 
the  integral  equations.  An  interpolation  between  the  results  of 
these  twe  methods  gives  very  good  results  over  the  entire  fre¬ 
quency  range,  even  when  few  basis  functions  are  used.  Results 
for  a  variety  of  configurations  are  shown  and  are  compared  with 
experimental  data  and  other  numerical  techniques. 


I.  Introduction 


which  was  extended  to  higher  frequencies  by  Kennelly  “ 
and  Affel  12].  Haefner’s  1937  paper  [3]  represents  the  ,-' 
most  extensive  experimental  data  on  the  resistance  of 
rectangular  conductors  with  a  wide  variety  of  width-to- 
thickness  ratios.  More  recently,  Weeks  et  al.  [4]  did  simi¬ 
lar  work  as  part  of  a  theoretical  treatment  of  the  problem. 

In  terms  of  theoretical  work,  the  circular  conductor  was  , 
the  first  case  considered,  since  it  allows  an  analytical 
solution  Maxwell  (5]  examined  nonperiodic  current; 
Kelvin  [6]  solved  the  periodic  case.  Carson  [7]  gave  a 
series  solution  for  the  two-wire  proximity  effect.  Cock¬ 
croft  I8)  used  the  Schwarz-Chrisloffel  transformation  to 
obtain  a  high-frequency  approximation  to  the  skin  effect 
which  was  expressed  in  terms  of  elliptic  integrals.  Wheeler  . 
[9]  discussed  the  “incremental  inductance”  rule,  which  is 
a  high-frequency  estimate  of  both  the  skin  and  proximity 
effects.  More  recently,  Casimir  and  Ubbink  [10],  (11) 
presented  an  overview  and  summary  of  the  basics  of  the 
skin  effect,  with  formulas  for  the  high-frequency  limits  of 
simple  cases. 

In  the  “filament  technique,”  the  conductor  (usually 
rectangular)  is  divided  into  a  large  number  of  rectangular 
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WITH  the  ever-increasing  speed  and  density  of  mod¬ 
em  integrated  circuits,  the  need  for  electromag¬ 
netic  wave  analysis  of  phenomena  such  a»  the  propagation 
of  transient  signals,  especially  the  distortion  of  signal 
pulses,  becomes  crucial.  One  of  the  most  important  causes 
of  pulse  distortion  is  the  frequency  dependence  of  con¬ 
ductor  loss,  which  can  be  incorporated  into  circuit  models 
for  transmission  lines  as  frequency-dependent  resistance 
and  inductance  per  unit  length.  Experimental  work  mea¬ 
suring  the  resistance  and  inductance  of  conductors  has 
concentrated  on  circular  and  rectangular  cross  sections. 
Kennelly  et  al.  [1]  did  a  thorough  experimental  study. 


filaments,  which  are  considered  to  have  uniform  current 
distribution  within  them.  Graneau  (12]  uses  a  power-series 
approach  in  frequency,  which  Weeks  el  al.  [4]  dispensed 
with.  Silvester  expands  the  current  in  a  flat  conductor  [13] 
in  a  series  of  eigenmodes  and  the  current  in  a  conductor 
of  arbitrary  shape  [14]  in  filaments.  In  both  cases  he 
ignores  the  effect  of  the  placement  of  the  return  current, 
or,  in  other  words,  the  proximity  effect.  While  these 
filament  methods  tend  to  be  very  good  for  low  frequen¬ 
cies,  since  the  current  density  is  then  almost  uniform,  they 
do  not  model  singularities  of  the  current  density  at  high 
frequencies  well. 
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The  Other  class  of  methods  involves  solving  the  mag¬ 
netic  vector  potential  integral  equation  [15)-i2l].  The 
boundary  condition  is  usually  on  the  tangential  magnetic 
field,  specified  on  a  closed  contour  some  distance  from 
the  conductor.  The  integral  equation  is  solved  by  tradi¬ 
tional  matrix  methods,  cither  by  expanding  the  current  in 
a  series  of  orthogonal  eigenfunctions  (15]-(18]  or  by  divid¬ 
ing  the  current  into  subdomain  basis  functions  (19]-(21]. 
This  method  is  limited  in  that  it  requires  knowledge  of 
the  magnetic  field  outside  the  conductor  somewhere  to 
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calculate  the  current  distribution  inside  but  does  not  give 
any  general  way  to  determine  that  field.  The  more  recent 
work  of  Cangellaris  [22]  applies  the  boundary  conditions 
developed  in  the  filament  approaches  to  the  magnetic 
vector  potential  integral  equation,  thus  removing  one  of 
the  principal  difficulties  of  that  method;  however,  it  still 
requires  modeling  of  the  current  throughout  the  cross 
section. 

In  recent  years,  new  methods  have  been  developed 
which  require  modeling  the  current  distribution  only  on 
the  surface  of  the  wires,  rather  than  throughout  the  cross 
section.  Djordjevic  et  al.  [23]  assumed  a  nonphysical  dis¬ 
tribution  of  current  along  the  propagation  direction,  which 
led  to  an  excess  resistance  at  high  frequencies.  Their  work 
was  modified  by  Wu  and  Yang  [24]  to  allow  appropriate 
quasi-TEM  propagation.  However,  since  both  of  these 
methods  depend  on  the  calculation  of  the  normal  deriva¬ 
tive  of  the  current  density,  they  have  numerical  difficul¬ 
ties  at  low  frequencies,  when  the  current  is  almost  uni¬ 
form  and  the  normal  derivative  is  small. 

The  technique  presented  in  this  paper  is  hybrid  cross- 
section  coupled  circuit/surface  integral  equation  ap¬ 
proach.  For  low  frequencies,  a  filament  method  based  on 
the  work  of  Weeks  et  al.  is  used,  except  with  triangular 
rather  than  rectangular  patches.  For  high  frequencies,  a 
surface  integral  equation  method  is  used.  However,  in 
contrast  to  previous  work,  the  calculation  of  resistance 
and  inductance  is  based  on  power  dissipation  and  stored 
magnetic  energy,  rather  than  on  impedance  ratios.  It  will 
therefore  be  more  easily  extended  to  structures  where 
nonuniform  propagation  can  occur.  In  the  middle  range 
of  frequency,  an  interpolation  is  made  between  the  results 
of  the  two  methods.  Since  this  is  a  frequency-domain 
method,  we  will  assume  an  e~“’'  dependence  to  all  quan¬ 
tities. 


II.  Cross-Section  Coupled  Circuit  Method 


distribution  is  used.  The  resistance  and  inductance  matri¬ 
ces  for  the  patches  (r  and  /)  are  then  calculated,  where 
these  matrices  are  defined  by 
di> 

—  =  (/a,/-^)i  (2) 

where  v  is  the  column  vector  of  the  voltage  differences 
between  the  patches  and  the  reference  patch,  and  i  is  the 
column  vector  of  currents  flowing  in  the  f  direction 
through  the  patches.  There  are  two  conditions  on  the 
system:  first,  that  the  total  current  in  each  wire  be  the 
sum  of  the  currents  in  the  patches  and,  second,  that  the 
voltage  on  each  patch  in  a  wire  be  the  same,  since  no 
transverse  currents  are  allowed  under  the  quasi-TEM 
assumption.  Using  these  conditions,  the  matrices  for  the 
patches  can  be  reduced  to  the  matrices  for  the  wires. 

For.  the  calculation  of  r  and  /,  we  follow  [4]  quite 
closely.  The  elements  of  the  resistance  matrix  of  the 
patches  are 

1  1 

1 

0*.mn  =  — 7-.  (3) 

where  the  first  subscript  indicates  the  wire,  the  second 
the  patch  within  the  wire;  is  the  cross-sectional  area 
of  patch  k  on  wire  j,  and  patch  0  on  wire  0  is  the 
reference  of  voltage.  Also  following  [4],  the  elements  of 
the  inductance  matrix  can  be  written  as  the  sum  of  partial 
inductances: 

~  ljk!mn  ~  ^>*.00  “  ^OO^mn  ^OO.'oO  (^) 

where  the  partial  inductances  are  given  by 


itp) 

*Jk,mn 


/  /  dS,,  / 1  dS'. 


For  low  frequencies,  we  use  a  two-dimensional  cross- 
section  coupled  circuit  method  to  find  the  resistance  and 
inductance  matrices  for  multiple  transmission  lines  with 
uniform  cross  sections.  We  assume  that  these  transmis¬ 
sion  lines  consist  of  signal  lines  over  a  common  return 
path  or  “ground  plane.”  The  matrices  R  and  L  are 
defined  by 

■  dV  =  ^ 

(') 

where  V  is  the  column  vector  of  the  voltage  differences 
between  the  wires  and  a  reference  wire  (ground  plane  or 
return  conductor),  and  /  is  the  column  vector  of  currents 
flowing  in  the  wires. 

Here  is  an  overview  of  the  cross-section  coupled  circuit 
method.  Each  conductor  is  divided  into  triangular  patches 
and  one  of  the  patches  from  the  return  conductor  is 
chosen  to  be  the  reference.  The  current  is  assumed  uni¬ 
form  on  the  cross  section  of  each  patch;  in  other  words,  a 
piccewise-constant  approximation  to  the  actual  current 


•ln[(x-xf+(y-y')^]  (5) 


where  x  and  y  are  coordinates  on  patch  jk,  and  x'  and 
y'  are  coordinates  on  patch  mn. 

In  the  Weeks  method,  the  patches  over  which  the 
integrals  in  (5)  are  done  are  rectangles,  and  the  quadruple 
integral  is  done  quite  easily  in  closed  form.  However,  if  is 
also  possible  to  evaluate  the  quadruple  integral  in  closed 
form  for  any  polygonal  shapes;  the  details  are  rather 
complex  and  are  left  for  the  Appendix.  We  therefore  use 
triangular  patches  as  the  most  flexible  means  of  modeling 
conductors  with  arbitrary  cross  sections;  polygons  are 
covered  exactly,  and  we  are  able  to  model  quite  closely 
other  shapes,  such  as  circles. 

Once  the  resistance  and  inductance  matrices  for  the 
patches  have  been  obtained,  we  proceed  in  the  following 
manner.  Taking 


dv 

dz 


r 


(6) 


The  matrix  ^  is  inverted,  and  =  Writing  out  the 
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elements  of  the  matrix, 

1;  in 

m -0  n -  1  ^ 

where  N  is  the  number  of  wires  and  N„  is  the  number  of 
patches  on  wire  m.  The  conditions  on  v  and  i  discussed 
above  are  applied,  to  give 

dV 

(8) 

m  -  1 

where  V  and  /  are  the  voltage  and  current  column  vectors 
for  the  wires  and  where 

E  Ey>*.«n-  (9) 


Inverting  y  givs 


Y-'^R-iwL.  (10) 

Thus,  the  frequency-dependent  resistance  and  inductance 
matrices  for  the  wires  have  been  obtained. 

In  [4],  the  distribution  of  patches  was  a  function  of 
frequency;  as  the  frequency  increased,  the  patches  were 
concentrated  at  the  edges,  where  the  current  is.  However, 
as  shall  be  shown,  it  is  more  efficient  to  switch  to  a 
surface  integral  equation  technique  for  high  frequencies; 
in  this  paper  the  distribution  of  triangular  patches  is  not 
altered  as  the  frequency  is  increased.  This  has  the  advan¬ 
tage  that,  since  the  resistance  and  inductance  matrices  of 
the  patches  are  independent  of  frequency,  r  and  I  need 
be  calculated  only  once,  no  matt^  for  how  many  frequen¬ 
cies  we  wish  to  calculate  R  and  L. 


Ill.  Differential  Equations  and 
Boundary  Conditions 

In  this  section,  the  basic  equations  and  boundary  condi¬ 
tions  which  will  be  used  in  the  surface-integral  equation 
method  will  be  derived.  The  coordinate  system  used  is 
shown  in  Fig.  1.  We  will  rely  heavily  on  quasi-TEM 
assumptions.  First,  outside  the  wires,  we  assume  that  the 
fields  are  transverse,  and  that  they  obey  Laplace’s  equa¬ 
tion.  In  other  words.  Maxwell’s  equations  outside  the 
wires  become 

Vxff  =  0  (11) 

Vi/  =  0  (12) 

where  V  here  is  only  the  transverse  operator,  xd/dx  + 
yd /dy.  The  vector  magnetic  potential.  A,  is  defined  such 
that  fiH  =  V  X  /4.  By  the  quasi-TEM  assumption,  J zJ^ 
and  A=^  zA^,  and  it  can  be  shown  that 

=  (13) 

Also,  inside  the  conductors,  the  displacement  current  is 
ignored;  Maxwell’s  equations  become 

VxE’^iwfiH  (1^) 

vx//  =  y  (15) 

V//  =  0.  (16) 


Fig.  1.  Coordinate  system  for  surface  integral  equation  method. 


Using  J  <*=  <tE,  this  reduces  to 

W, -I- iwAio-y,  ■=  0.  (17) 

There  are  two  boundary  conditions  at  the  interface  of 
the  conductors  and  free  space:  the  continuity  of  tangen¬ 
tial  H  and  of  normal  B  =  fiH.  If  H  outside  is  expressed 
in  terms  of  A and  H  inside  in  terms  of  J^,  the  condition 
on  H  reduces  to 


which  is  satisfied  along  all  the  conductor-free  space  in¬ 
terfaces.  The  boundary  condition  on  normal  B  is  more 
difficult.  Assuming  that  all  the  materials  have  the  same 
permeability. 


If  the  derivatives  of  two  quantities  along  a  line  are  equal, 
then  those  quantities  must  be  equal,  to  within  a  constant: 

y, - /!,].  (20) 

is  constant  over  a  single  conductor,  but  can  vary  from 
conductor  to  conductor. 

Finally,  it  is  necessary  to  be  able  to  specify  the  total 
current  flowing  on  a  wire.  Using  Green’s  first  identity, 

/|d5(d>VV  +  Vd.-Vi^)=^d/d)^  (21) 

where  dS  is  a  two-dimensional  integral  over  a  cross-sec¬ 
tional  area,  and  dl  is  a  one-dimensional  integral  along  the 
closed  contour  bounding  that  area.  Also,  the  normals  are 
defined  as  pointing  out  from  the  region  of  interest.  With 
lA  “  y,  and  ”  1,  and  considering  (17), 

I~  f  fdSJ,~  —  f  fdSV^J^  =  —6dl^  (22) 
which  is  an  expression  for  the  total  current  flowing  in  a 
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wire  in  the  z  direction  in  terms  of  quantities  on  the 
surface. 

IV.  Determination  of  Circuit  Parameters 

Eventually,  the  quantities  of  interest  are  the  resistance 
and  inductance  per  unit  length  of  these  conductors.  It 
turns  out  that  it  is  possible  to  express  these  quantities  in 
terms  of  the  current  and  its  normal  derivative  on  the 
surface  of  the  wires.  This  is  useful,  because  it  means  that 
the  problem  can  be  formulated  in  terms  of  a  set  of 
coupled  integral  equations  involving  only  these  surface 
quantities,  allowing  a  great  savings  in  computation.  The 
resistances  and  inductances  will  be  derived  through  power 
and  energy  considerations. 

For  the  resistance,  consider  a  case  with  a  current  I 
flowing  in  a  signal  wire  and  returning  in  a  reference,  for 
example  a  ground  plane.  Starting  with  the  power  defini¬ 
tion  of  resistance, 

2P  fjdSE.J:  j  jjdS\jf 

R - ^  =  — - =  -  — -  (23) 

>^1  \jJdSJf 

where  the  integration  in  the  numerator  is  over  all  wires, 
while  that  in  the  denominator  is  only  over  the  signal  wire. 
The  numerator  can  be  put  in  a  more  useful  form,  using 
(21)  with  <f»  =  y,  and  i{i  -=  J*,  and  its  complex  conjugate, 
together  with  (17)  to  get 


side  the  wires. 


ffds\jf~f{dsjj: 


2  u/jicr 


1  /  ,  r  dj* 
- — 

2  [  cn  an 


Including  the  total  current  squared,  |/1^,  from  (22), 


R  =  ii>ix 


6  dllmlj,  —  ] 

~at\  wires _ (  dn  ) 

r  ' 

6  di-^ 

''signal  wire 


X  ( 

W  dl  RcudtcrA^— — 

^11  wires _ I _ 

6  di-^ 

‘'signal  wire 


The  mutual  resistance  and  inductance  are  calculated 
from  energy  considerations,  and  from  the  self  terms  calcu¬ 
lated  above.  If  we  specify  a  current  to  flow  on  line  /, 
and  —  to  flow  on  line  j,  we  can  calculate  the  power 
dissipated,  P^,  and  the  stored  magnetic  energy,  W„,  very 
easily  by  the  above  technique.  This  gives 

=  (28) 


(29) 

V.  Derivation  of  Coupled  Integral  Eouations 

In  order  to  complete  the  formulation,  integral  equa¬ 
tions  are  required  which  relate  A^to  dA^  /dn  outside  the 
wires  and  to  dj^/dn  inside  the  wires.  Starting  from 
Green’s  theorem: 

,  ,  /.(  dip  dd>  \ 

j /dS'(.^VV-'AV'<f»)=^d/j<^— -(A— J  (30) 

where  the  integral  dS'  is  over  a  cross-sectional  region,  the 
integral  dV  is  over  the  contour  bounding  that  region,  and 
the  normals  point  out  from  the  region  of  interest.  In 
general,  let  'I'(p)  be  either  A^  or  where  p  is  position 
in  the  two-dimensional  cross  section.  Since  ^(p)  satisfies 
+  C’P  =  0,  where  C  =  0  for  Laplace’s  equation  and 
C  =  for  the  diffusion  equation,  a  Green’s  function, 
G(p,  p’),  can  be  found  which  satisfies  V^G -t- CG  = 
-  5(p  -  p').  Substituting  'P  for  ip  and  G  for  <f)  in  (30), 

jldS'np')Hp-p) 

f  ^'A’(p')  ^G{p,p') 

--i^dr  G{p,p')-^-^{p')—-  .  (31) 


Similarly,  starting  with  a  magnetic  stored  energy  defini¬ 
tion  of  inductance  per  unit  length,  L: 

UdSHH* 

- - -  (26) 

1/1'  1/1' 

where  the  range  of  integration  for  the  numerator  is  over 
all  space.  Using  a  technique  similar  to  that  for  (25), 
combining  the  contributions  from  regions  inside  and  out¬ 


The  integral  equation  will  be  formulated  on  the  surface, 
so  both  p  and  p'  are  on  the  surface.  This  places  6(p  -  p’) 
just  on  the  boundary  of  the  dS'  integration,  and  this  must 
be  treated  carefully.  The  most  straightforward  method  is 
that  integrating  a  delta  function  that  lies  on  the  edge  of 
the  range  of  integration  gives  1/2.  With  this. 


^^di'nn  + 

Integral  equations  for  both  A,  and  J.  can  now  be 
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obtained.  For  the  outer  equation  (13), 


(f) 

■'all  wires 


-6  di'A.in 

‘'•II  wires 


\dG„il,r)  1 

■)  -jHi-n  (33) 


where 


G„(p.p')--:^ln 


v/(ar-x')'  +  (y-yf  ]  (34) 


and  where  the  sign  change  in  (33)  is  due  to  the  normals 
pointing  into  the  region  of  interest.  The  range  of  integra¬ 
tion  is  over  the  boundary  of  the  free-space  region,  in 
other  words,  over  the  surface  of  every  wire. 

Similarly,  for  the  inner  equation  (17),  for  each  wire. 


•'wireo  Oft 


j.  r^G,(/,/')  1 

d/'7,(/')  — ^  +  -5(/-/')  (35) 

^wire^  [  oft  Z 


where 


where  “ker”  and  “kei"  are  the  real  and  imaginary  Kelvin 
functions.  This  equation  applies  to  each  wire  separately; 
the  range  of  integration  is  over  the  surface  of  that  wire. 

Using  the  boundary  conditions  (18)  and  (20)  to  elimi¬ 
nate  A,  from  the  integral  equation  for  the  outside  fields 
and  add  the  condition  on  the  total  currents  from  (22), 

(f>  dVGSi,n-^-<b  dr 

■'jU  wires  otl  'all  wires 

•(y,(/')-hiW/i,]  -^^-r-^--g(/-/')  (37) 


There  is  one  important  thing  to  note  about  these  inte¬ 
gral  equations.  Contrary  to  simpler  electrostatic  prob¬ 
lems,  the  boundary  conditions  are  neither  Dirichlet  nor 
Neumann;  both  G  and  dG/dn’  must  be  kept  in  the 
equations.  Because  of  this,  the  formulation  is  in  terms  of 
the  free-space  Green’s  functions  for  Laplace’s  equation 


and  the  diffusion  equation,  (34)  and  (36).  All  information 
about  the  boundaries  is  contained  in  the  paths  of  integra¬ 
tion. 


VI.  Solution  of  Coupled  Integral  Equations 

We  solve  these  coupled  integral  equations,  (37),  (38), 
and  (39),  by  the  method  of  moments  with  subdomain  basis 
functions.  Expanding  the  unknown  functions  7,  and 
SJj/dn  as  the  sum  of  known  functions  times  unknown 
coefficients, 

(40) 

m 

^“E*„B„(/).  (41) 

nt 

Simple  pulse  basis  functions  are  used,  normalized  so  that 
the  integral  is  unity: 


^  n/A„  if/„</</„-FA 
\o  otherwise. 


(42) 


This  gives  a  piecewise-constant  approximation  to  the  sur¬ 
face  quantities.  The  same  functions  are  used  for  testing, 
thereby  implementing  Galerkin’s  method. 

It  turns  out  that,  for  high  frequencies,  the  current 
distribution  on  a  wire  is  similar  to  the  charge  distribution 
on  a  perfect  conductor.  For  polygonal  wires,  this  means 
that  the  current  will  be  concentrated  at  the  comers. 
Theiefore,  we  find  it  advantageous  to  concentrate  the 
basis  functions  in  the  same  way.  For  three  basis  functions 
on  a  side,  for  example,  the  two  in  the  corners  are  each 
one  eight  the  length  of  the  side;  the  center  one,  three 
quarters.  These  values  were  determined  empirically,  by 
seeing  which  division  gave  results  closest  to  those  for  a 
large  number  of  basis  functions. 

We  can  thus  approximate  the  coupled  integral  equation 
as  a  matrix  equation: 


where  J  is  the  vector  of  the  unknown  /„’s  (current),  K  is 
the  vector  of  the  unknown  (normal  derivative  of  the 
current),  and  Ag  is  the  vector  of  the  A^'%  (constants  of 
vector  potential).  The  total  currents  on  e^h  wir^  a£e 
specified  by  the  vector  I.  The  matrices  Vg,  iVg,  Ug,  J,  F], 
and  Uf  arise  from  integrals  of  products  of  the  Green’s 
functions  with  the  basis  functions  and  are  completely 
known.  The  solution  of  this  matrix  equation  by  LU  de¬ 
composition  provides  us  with  an  approximation  for  J,  and 
dj,/dn,  and  through_(25)  and  (27^  R  and  L.  Since  the 
outer  matrices  ,  Ug,  IVg,  and  S,  are  independent  of 
frequency,  they  only  need  to  be  calculated  once.  We  can 
make  use  of  this  fact  by  LU-decomposing  this  part  of  the 
large  matrix  only  once.  Completing  the  decomposition 
with  the  rest  of  the  matrix  for  each  value  of  frequency. 
For  the  pulse  basis  functions  used,  the  outer  matrices  can 
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Fig.  2.  Resistance  of  isolated  square  wire. 


^  expressed  in  closed  form.  The  inner  matrices,  and 
Ui,  have,  at  worst,  double  numerical  integrals,  and  for  the 
case  of  the  interaction  between  elements  which  lie  on  the 
same  line,  these  can  be  expressed  in  closed  form.  Also,  at 
high  frequencies,  owing  to  the  highly  local  nature  of  the 
diffusion  Green’s  function  (36),  which  results  from  the 
rapidly  decaying  asymptotic  nature  of  the  Kelvin  func¬ 
tions,  only  the  neighboring  patches  have  an  appreciable 
interaction;  those  integrals  can  be  calculated  quite  rapidly. 

VII.  Results 

In  these  results,  we  will  compare  the  hybrid  method, 
described  above,  with  experimental  results,  as  well  as  with 
two  other  methods:  the  Weeks  method  [4],  which  models 
the  current  throughout  the  cross  section,  and  the  work  of 
Djordjevic  et  al.  (23),  which  models  an  equivalent  current 
only  on  the  surface  over  all  frequencies.  As  shall  be  seen, 
by  using  a  hybrid  method,  we  can  avoid  the  weaknesses  of 
both  of  these  methods. 

First,  we  consider  the  example  of  an  isolated  square 
conductor,  4.62  mm  on  a  side,  with  conductivity  <r  —  5.72 
X  10’  (ft  -m)"'.  While  the  inductance  per  unit  length  is 
undefined,  the  hybrid  method  can  be  used  to  calculate 
the  resistance  per  unit  length  and  for  comparison  with  the 
experimental  results  of  Haefner  [3]  and  the  results  ob¬ 
tained  by  using  the  Weeks  method  [4],  As  can  be'  seen 
(Fig.  2),  the  fit  for  the  new  method  is  quite  good.  In  this 
example,  only  12  basis  functions  (25  unknowns)  were 
used,  three  on  a  side.  By  comparison,  the  Weeks  method 
with  49  basis  functions  does  not  give  as  good  a  result. 

The  next  example  is  that  of  two  parallel  circular  wires, 
a -5.84x10’  (fl-m)"'.  The  wires  have  a  diameter  of 
11.68  mm,  and  a  separation  of  0.3  mm  and  8  mm  for  the 
two  cases.  Here,  the  circles  were  modeled  as  n-sided 
polygons  having  the  same  cross-sectional  area.  In  Figs.  3 


Frequency  (Hz) 

Fig.  3.  Resistance  of  two  circular  wires. 


and  4,  the  results  for  13  basis  functions  per  circle  (54 
unknowns)  are  shown,  compared  with  the  experimental 
results  of  Kennelly  et  al.  [1],  The  fit  is  again  quite  good 
with  the  experimental  results.  Since  the  Weeks  method  is 
limited  to  rectangular  elements,  it  is  not  capable  of  han¬ 
dling  this  case. 

Next,  we  take  the  example  of  two  parallel  rectangular 
wires,  o-  —  5.6  x  10’  (.Cl  -m)" the  configuration  is  shown 
in  Fig.  5.  In  Figs.  6  and  7,  we  compare  the  hybrid  method, 
the  Weeks  method,  and  the  results  from  [23]  calculated 
from  a  purely  surface  integral  equation  approach.  It  can 
be  seen  that  the  hybrid  method  agrees  with  each  of  the 
others  in  its  range  of  validity.  Also,  the  numerical  instabil¬ 
ity  of  purely  surface-integral  equation  methods  in  calcu- 
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2  mm  2  mm 


Case  I  Case  il 

Fig.  S.  Two  rectangular  wires. 


Frequency  (Hz) 


Fig.  6.  Resistance  of  two  rectangular  wires. 

lating  the  low-frequency  inductance  can  be  observed.  It  is 
also  clear  that  the  hybrid  method  in  general  requires 
fewer  basis  functions,  and  thus  less  computation  time, 
than  the  Weeks  method.  In  fact,  as  the  frequency  in¬ 
creases  and  the  conductors  become  many  skin  depths 
across,  even  a  large  number  of  basis  functions  in  the 
Weeks  method  leaves  us  with  a  significant  error.  This  is 
due  to  the  inability  of  the  Weeks  method  to  correctly 
model  the  distribution  of  current  along  the  surface,  which 
is  crucial  to  the  calculation  of  resistance  at  such  frequen¬ 
cies.  Tables  1  and  II  compare  the  results  of  the  hybrid 
method  with  the  Weeks  method  for  the  case  of  two 
square  wires,  including  CPU  times,  on  a  Digital  Equip¬ 
ment  Corporation  VAXstation  3500,  running  VMS.  As 
can  be  seen,  the  cost  of  the  hybrid  method  in  terms  of 
CPU  time  is  much  lower  than  the  Weeks  method  for 
anything  more  than  a  moderate  number  of  basis  func¬ 
tions,  and  especially  for  high  frequencies. 

Finally,  we  consider  the  case  of  three  rectangular  con¬ 
ductors  over  a  ground  plane,  <7  •=  5.81x10’  (fl— m)~'. 
The  configuration  is  shown  in  Fig.  8,  the  resistance  of  the 
first  line  in  Fig.  9,  and  the  self-  and  mutual  inductances  in 
Fig.  10.  Since  our  method  is  not  capable  of  modeling  an 
infinite  ground  plane  (since  imaging  is  very  difficult  with 
an  imperfect  conductor),  a  very  large  conductor  was  used, 
4  mm  by  0.5  mm,  in  its  place,  with  more  basis  functions 
concentrated  in  the  area  under  the  signal  lines.  We  com¬ 
pare  with  the  Weeks  method  and  with  a  purely  surface 


Frequency  (Hz) 


Fig.  7.  Inductance  of  two  rectangular  wires. 


Fig.  8.  Three  rectangular  wires  over  a  ground  plane. 


integral  equation  result,  both  with  a  large  number  of  basis 
functions.  The  same  trend  can  be  seen  as  in  the  previous 
case,  but  there  is  now  a  large  error  in  the  high-frequency 
inductance  as  predicted  by  the  Weeks  method.  This  is 
due  to  the  fact  that  the  Weeks  method  does  not  model 
the  concentration  of  the  current  on  the  ground  plane 
under  the  signal  lines,  since  at  high  frequencies  most  of 
the  patches  are  concentrated  at  the  corners  of  the  ground 
plane,  far  away  from  the  current.  If  one  improves  the 
Weeks  method  by  restricting  the  majority  of  basis  func¬ 
tions  to  be  under  the.signal  lines,  one  gets  results  which 
agree  with  the  hybrid  method  quite  closely. 

VIII.  Conclusions 

A  technique  has  been  developed  to  calculate  the  skin 
effect  resistance  and  inductance  of  transmission  lines  with 
arbitrary  cross  sections.  This  technique  provides  accurate 
answers  over  a  wide  range  of  frequencies,  including  the 
range  where  neither  low-frequency  (direct  current,  uni¬ 
form  distribution)  nor  high-frequency  (skin  depth)  ap¬ 
proximations  are  valid.  The  technique  is  a  hybridization 
of  two  distinct  methods.  The  first  is  a  cross-section  cou¬ 
pled  circuit  approach,  subdividing  the  wires  into  triangu¬ 
lar  patches  which  are  assumed  to  have  uniform  current 
distribution.  This  method  is  best  for  low  frequencies, 
when  the  physical  current  has  very  little  variation  across 
the  cross  section.  The  second  method  is  in  terms  of  a 
coupled  integral  equation,  linking  the  current  and  its 
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TABLE  1 

Results  and  CPU  Times  for  Two  Square  Wires  with 
Hybrid  Method 


Basis 

Functions 

Number  of 
Unknowns 

Frequency  R  L 

(Hz)  (mfl/m)  (nH/m) 

CPU  Time 
(s) 

3x3 

50 

(Preprocessing) 

12.99 

10’ 

8.929 

599.5 

1.35 

10* 

11.07 

577.7 

3.36 

10‘ 

94.76 

466.6 

1.24 

7x7 

114 

(Preprocessing) 

19.16 

10’ 

8.929 

599.5 

1.37 

10* 

11.13 

579.6 

17.00 

10‘ 

98.54 

466.5 

6.77 

ISxlS 

242 

(Preprocessing) 

155.79 

10’ 

8.929 

599.5 

33.44 

10* 

11.15 

580.2 

156.49 

10* 

98.84 

466.9 

50.06 

TABLE  II 

Results  and  CPU  Times  for  Two  Square  Wires  with 
THE  Weeks  Method 


Basis 

Functions 

Number  of 
Unknovms 

Frequency 

(Hz) 

R 

(mfl/m) 

L 

(nH/m) 

CPU  Time 
(s) 

3x3 

17 

10’ 

8.929 

599.5 

0.76 

10* 

10.46 

586.1 

0.76 

10* 

118.8 

466.9 

0.75 

7x7 

97 

10’ 

8.929 

599.5 

50.14 

10* 

11.10 

581.1 

49.87 

10* 

92.60 

468.5 

49.43 

15x15 

449 

10’ 

8.929 

599 J5 

3684.42 

10* 

11.22 

580.1 

3630.31 

10* 

91.80 

468.1 

3639.78 

100 


E 

S  10 


-  hybrid  Technique  (104) 

- Surfote  Int.  Tech.  (536) 

O  Weeks'  Method  (538) 

+  Improved  Weeks  (167) 


1000 


10^ 


10’ 


10® 


10' 


10® 


10” 


Frequency  (Hz) 


Fig.  9.  Resistance  of  ihree  rectangular  wires  over  a  ground  plane. 


Fig.  10.  Self-  and  mutual  inductances  of  three  rectangular  wires  over  a 
ground  plane. 


normal  derivative  on  the  surface  of  each  wire  with  the 
magnetic  vector  potential  and  its  normal  derivative  on  the 
same  surfaces;  the  resistance  and  inductance  are  both 
expressed  in  terms  of  these  surface  quantities.  This 
method  is  best  for  high  frequencies,  when  the  current  is 
almost  all  confined  to  the  surface,  and  the  diffusion 
Green’s  function  (eg.  (36))  is  very  localized.  For  the 


middle  frequency  range,  an  interpolation  between  the  two 
results  gives  very  good  accuracy  with  few  basis  functions. 
The  interpolation  function  was  based  on  the  average  size 
of  the  conductors,  measured  in  skin  depths,  and  was  of 
the  form  1/(1  +  0.16fl^/5^),  where  a  is  the  average  cross 
section  of  the  conductors,  and  5  is  the  skin  depth.  The 
optimization  of  the  interpolation  function  is  an  area  of 
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further  research.  By  choosing  triangular  patches  for  the 
cross-section  method  and  free-space  Green’s  functions  for 
the  surface  method,  a  single  program  is  able  to  handle 
arbitrary  conductors.  The  method  is  limited  at  present  to 
infinite,  uniform  lines,  although  nothing  theoretically  pro¬ 
hibits  extension  to  three-dimensional  lines.  The  theory 
behind  this  method  is  not  nece.ssarily  limited  to  configura¬ 
tions  with  uniform  dielectrics,  but  problems  in  the  de¬ 
finitions  of  resistance  and  inductance,  stemming  from 
difficulties  with  the  extension  of  current  and  especially 
voltage  to  non-TEM  lines,  make  such  an  extension  of  the 
method  not  immediately  obvious.  For  most  practical  cases, 
however,  the  effects  of  nonuniform  dielectrics  on  the 
resistance  and  inductance  can  be  ignored,  so  that  the 
method  presented  in  this  paper  will  give  quick  and  accu¬ 
rate  results. 


Appendix 

Closed-Form  Expression  for  Partial  Inductances 
The  problem  is  to  evaluate  the  integral 


I-jjdS„jjdS‘,„\nt  (Al) 


where  f  e(x  -  x')^  +(y  -  y')^  and  the  areas  of  integra¬ 
tion  are  the  triangular  patches  (y)  and  (km).  Using  the 
fact  that  In r  =  V^V'^t^dn t -3)/64,  and  Green's  first 
identity  (21) 


where  the  integrals  dl,^  and  dl't.„  are  over  the  perimeters 
of  patches  (ij)  and  (km),  respectively,  and  the  normals 
point  out  from  the  patches.  Using  the  chain  rule. 


I  =  ■^j>dl^,<f)dll 


ah 

- 1 

inf-- 

dn  dn' 

2 

■¥ 


dt  dt  (  3 

- r  In'-T 

dn  dn  \  2 


(A3) 


If  the  patches  are  polygons,  these  integrals  over  the 
perimeters  of  the  wires  become  just  sums  of  integrals  over 
pairs  of  line  segments.  Therefore,  we  need  to  be  able  to 
evaluate  this  integral  where  the  paths  of  integration  are 
arbitrarily  oriented  line  segments.  Without  loss  of  gener¬ 
ality,  the  coordinate  system  is  chosen  so  that  the  dl 
segment  is  parallel  to  the  x  axis  (Fig.  11).  In  this  case,  the 
various  derivatives  of  t  are 


r»(u-Fx„-  X, cos  +  ( y^  -  y,,- sin  <f>v)^  (A4) 

at  at 

— — =  2(y,-y„-t-sin0u)  (A5) 

3n  ay,, 

at  dt  at 

—  “  -  cos  <t> - -(■  sin  <fii - 

an  ay,  dx,, 

“  -2[sind>(u  -  X,  +  xJ+cos<f>(y,.-y„)  (A6) 


Fig.  II..  Coordinate  system  for  mutual  inductance  of  triangular  patches. 


and 


a^t  d  (  at  \ 

andn'  “  ~ay^^[a^j  “ 


(A7) 


Letting  x  s  x^.  -  x„  and  ysy,,  —  y„,  the  following  double 
integral  over  the  pair  of  line  segments  is  obtained: 


f(u,u)  j dvt (5  -l\n  t)  cos  4> 


-h(6-41n/)(y  -I-  Dsint^) 

•((u  -  x)sin</> -f- y  cos<^)  (A8) 


where  /  >=  (u  —  x -cost^u)^ -»-(y -fsin  ^t;)^.  These  inte¬ 
grals  can  be  done  in  closed  form  [25].  If  the  length  of  the 
dl  segment  is  a,  and  that  of  the  dl'  segment  is  b,  the 
contribution  to  (Al)  from  this  pair  of  line  segmenu  is 
/(fl,6)- /(fl,0)- /C0,f>)-i- /(O.O).  TTie  total  is  thus  the 
sum  of  the  contributions  from  each  pair  of  line  segments, 
one  from  the  (ij)  patch  and  the  other  from  the  (km) 
patch. 
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